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ABSTRACT 

The effective use of active magnetic bearings for vibration control in 
turbomachinery depends on an understanding of the forces available from a magnetic 
bearing actuator. The purpose of this project was to characterize the forces as functions 
shaft position. 

Both numerical and experimental studies were done to determine the characteristics 
of the forces exerted on a stationary shaft by a magnetic bearing actuator. The numerical 
studies were based on finite element computations and included both linear and nonlinear 
magnetization functions. 

Measurements of the force versus position of a nonrotating shaft were made using 
two separate measurement rigs, one based on strain guage measurement of forces, the 
other based on deflections of a calibrated beam. 

The general trends of the measured principal forces agree with the predictions of the 
theory while the magnitudes of forces are somewhat smaller than those predicted. Other 
aspects of theory are not confirmed by the measurements. The measured forces in the 
normal direction are larger than those predicted by theory when the rotor has a normal 
eccentricity. 

Over the ranges of position examined, the data indicate an approximately linear 
relationship between the normal eccentricity of the shaft and the ratio of normal to principal 
force. The constant of proportionality seems to be larger at lower currents, but for all cases 
examined its value is between 0. 14 and 0. 17. The nonlinear theory predicts the existence 
of normal forces, but has not predicted such a large constant of proportionality for the 
ratio. 

The type of coupling illustrated by these measurements would not tend to cause 
whirl, because the coupling coefficients have the same sign, unlike the case of a fluid film 
bearing, where the normal stiffness coefficients often have opposite signs. They might, 
however, tend to cause other self-excited behavior. This possibility must be considered 
when designing magnetic bearings for flexible rotor applications, such as gas turbines and 
other turbomachinery. 

In related work attached as an appendix, simulations of 2DOF systems subject to 
these force models show that significant nonlinear behavior can occur, including multiple 
coexisting solutions, bifurcations in response as the stabilities of the respective solutions 
change, and self-similarity in stability boundaries. 



11 


CONTENTS 

1 . INTRODUCTION 1 

1 . 1 Magnetic Suspension of Turbomachine Shafts 2 

1.1.1 Disadvantages 4 

1.1.2 Background 5 

1.1.3 Developments During this Project 6 

2 . TECHNICAL FINDINGS: Analytical/Numerical Modeling 9 

2.1 Fundamental Principle 9 

2.2 Air Gap Method 10 

2.3 Air Gap Method: Computer Program 14 

2.4 Air Gap Method: Results 15 

2.4.1 Test Case 1: parallel surfaces 15 

2.4.2 Test Case 2: Effects of element size and fringing 15 

2.4.3 Forces from one magnet of a bearing 15 

2.5 Full Magnet Method 18 

2.5.1 Differential Equations in 3-D 21 

2.5.2 Two-dimensional Equations 21 

2.5.3 Permeability 23 

2.5.4 Boundary Conditions 23 

2.5.5 Discretization 27 

2.6 Results of Linear Calculations 30 

2.7 Nonlinear Force Calculation 31 

2.7.1 Modelling of Magnetization Curve 31 

2.7.2 Calculation of Flux Distribution 32 

2.7.3 Results of Calculations 34 

2.8 Effects of Uncertainties and Property Variations 50 

2.9 Conclusions (Analytical/Numerical) 55 

3. TECHNICAL FINDINGS: Experiments 5 6 

3.1 Magnet Apparatus I 5 6 

3.1.1 Method of Measurement 60 

3.2 Results of Measurements, Apparatus 1 61 

3.3 Measurement Apparatus II 73 

3.3.1 Deflecting Beam Apparatus 73 

3.3.2 Measurement Method 75 

3.3.3 Measurements Using One Magnet 75 

3.3.4 Measurements Using Two Magnets 92 

4. CLOSURE 105 

4.1 Summary of Technical Findings 105 

4.2 Documentation 106 



1U 


CONTENTS, continued 

5 . REFERENCES 107 

APPENDICES (individually page numbered) 113 

APPENDIX A - PROGRAM FOR FORCE CALCULATION USING 

AIR GAPS ONLY 

APPENDIX B - COMPUTATIONAL METHODS INCLUDING METAL 

REGIONS 

APPENDIX C - MASTER’S THESIS OF THOMAS WALSH 



1. INTRODUCTION 

This report describes work done under Grant NAG 3-968 during the performance 
period October 1989 to 1 February 1992, in addition to further related work using 
knowledge gained during the work performed under this grant. The purpose of the 
research was to examine certain aspects of the potential of magnetic bearings for vibration 
control in turbomachinery. The principal thrusts of the research have been 

1) Calculation of the two-dimensional forces exerted on a shaft by a typical magnetic 

actuator under open loop conditions. 

2) Measurement of such forces using specially designed apparatus. 

3) Simulation of dynamics of a simple rotor using the measured and calculated forces 

along with a control law. 

The report consists of three principal sections, plus three appendices. Section 1 is 
an introduction and a brief review of pertinent literature at the time of the beginning of this 
work. 

Section 2 describes the analytical and numerical modelling of the magnetic flux 
distribution in magnets of a magnetic actuator and the results of calculations of force 
between the actuator and the shaft. Section 3 describes two kinds of experiments 
conducted to determine the forces that are modelled in Section 2. Because comparisons 
between theory and experiment are made, it is sometimes necessary to refer in Section 2 to 
force measurements that will be described in more detail in Section 3. Similarly, Section 3 
refers back to Section 2. 

Appendix A is a listing of the computer program used for calculations using an air- 
gap method. Appendix B is a description of methods used in calculations of force 
including flux contained in metal parts. 

Appendix C contains the text of a thesis submitted to Duke University by Thomas 
Walsh for the M.S. degree, which uses the results of measurements and calculations done 
under this grant in the simulation of a rotor-magnetic bearing system. This work was 
performed subsequent to the actual grant period, but is included because of the close 
relation to and dependence on the results obtained under this grant. 

In addition to the sections describing technical findings, this report summarizes 
related activities, including papers and reports, personnel, equipment and progress of 
students. 
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1.1 Magnetic Suspension of Turbomachine Shafts 


The introduction of practical magnetic supports for rotating shafts is a recent 
development. These devices have the potential of replacing fluid film bearings and rolling 
element bearings in some critical applications, and of acting as supplemental control 
actuators alongside these traditional bearings in order to limit vibration, noise and 
instability. Figure 1 . 1 shows the general concept of a rotating shaft suspended by a set of 
controlled magnetic actuators. In this report, the combination of an actuator, its sensors, 
controller and power amplifier is called a magnetic bearing. 

In a magnetic bearing the shaft is supported by the force established between a set of 
electromagnets and the shaft due to the magnetic field. There is no direct contact between 
the shaft and any part of the bearing. This method of support has several advantages over 
traditional fluid film or rolling element bearings. Since no lubricant is needed, there is no 
sealing requirement to prevent either the lubricant or the working fluid from contaminating 
the other. Also, the lubricant supply system required for a traditional bearing is eliminated, 
and the frictional losses in the magnetic support are negligible compared to those in a fluid 
film bearing. Finally, since the magnetic bearing requires a feedback control system to 
maintain stability even in a nonrotating steady-state case, this feedback loop may be used to 
advantage in adjusting the dynamic characteristics of the rotor-bearing system to optimize 
the machine s vibration characteristics. The magnetic bearing could be designed so that it 
opposes the destabilizing effects of other parts of the system. 

It is this last possibility that is the most exciting aspect of magnetic suspension. The 
ability to control better the dynamics of shafts and thus to reduce the danger and expense 
that result from high levels of vibration is the principal motivation for research in this field. 

Magnetic bearings are rapidly gaining acceptance as replacements for traditional fluid 
film bearings in the design of turbomachinery. Applications include the small and 
sensitive, such as turbomolecular pumps and x-ray generation equipment, as well as the 
extremely precise, such as machine tool spindles. At the other extreme, applications 
include very large industrial machinery such as compressors, turbines and engines [1], 

The motivations for these applications vary, but most are inspired by the possibility of 
precise control of the rotor through the magnetic bearing's active feedback loop. In the 
case of high precision machine tools there is an obvious need for precise control of tool 
position, and this control is made possible through the active magnetic bearing to a degree 
not possible by using other bearings, even the stiffest of rolling element bearings. This 
degree of control is possible even though the parameters of the bearing may not have been 
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Figure 1.1. The active magnetic bearing concept. 
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accurately assessed, because the active controller compensates for poorly known 
characteristics by the exercise of negative feedback. Thus, even though the system is not 
fully optimized, its performance may still far exceed that of passive bearings in terms of 
accuracy of positioning. 

The need for precision force generation is less obvious but no less crucial in the case 
of large industrial machinery. In this case, the machines are so large and expensive, and 
down-time is so costly, that it is essential to minimize vibration problems. Many of these 
machines operate at speeds higher than one or two critical speeds, so vibration problems 
can be intense, and the nature of these problems becomes successively worse as 
performance demands are increased. The traditional approach toward minimizing vibration 
is to design passive fluid bearings by choosing clearances and length-to-diameter ratios to 
achieve the best possible effective stiffness and damping characteristics. Most of the 
damping in these large vibrational systems arises from the bearings, and in passive 
bearings there is always a trade off between damping and stiffness. In addition, most such 
bearings give rise not only to principal restoring forces, generally desirable, but also to 
forces that act normal to a perturbation direction. Depending on their signs, these forces 
can be stabilizing or destabilizing. It is the nature of fluid bearings that in most practical 
situations they are destabilizing. Thus the active magnetic bearing, which offers 
controllable forces that are in theory uncoupled, has a strong appeal to the manufacturers 
and users of such machinery. 

1.1.1 Disadvantages 

Magnetic bearings are not a panacea, however. There are significant drawbacks to 
their use: magnetic bearings in general are larger and heavier than equivalent fluid 
bearings, they require continuous control and an uninterrupted power supply and therefore 
need redundant controllers and backup power capability, plus emergency backup bearings 
for shut-down in case of complete failure of the active system. These disadvantages must 
be weighed against the positive factors of reduced power losses, elimination of lubricant 
supply system and controllability. 

The disadvantages stemming from the size and weight factors may be minimized 
with better prediction of the forces available from actuators, and the effects of geometry on 
the available forces. Reliance on simplified theory for the forces available in an actuator 
has probably led to overdesign of the components, and to a reliance on controller 
robustness to compensate for inaccuracies in the force prediction. In practice the 
installation of a magnetic bearing system in a large turbomachine has been found to require 
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a lengthy process of tuning both on the test stand and later in the field for each individual 
machine 

Better force prediction will allow optimization of the actuators themselves in an open 
loop sense, shortening the tuning process and freeing the designer of the controls to 
concentrate on higher orders of vibration control strategy To this end, the work described 
in this report is concentrated on developing reliable and efficient methods of predicting the 
forces exerted by magnets on a rotor. The work consists of both theoretical/numerical 
analysis and experimental measurements of forces. 

1.1.2 Background 

This section describes the background and state of technology in magnetic bearings 
largely as it existed at the start of this project. The following section contains limited 
references to developments that occurred as this work proceeded. 

The concept of suspending a machine part by force of magnetic attraction was 
introduced as early as 1842 [2], and some early devices for magnetic support were 
attempted using permanent magnets and electromagnets, but practical application of the idea 
awaited relatively recent developments in control technology and power electronics. 

Beams [3] in 1949 built a successful magnetic suspension device for a small diameter 
rotor (1/64 inch) in order to achieve high rotational speeds. The system used vacuum tubes 
for control and power amplification, and thus was limited to supporting only small masses. 

The first application of fully active magnetic suspension was in the field of 
aerodynamic research where a system was developed to support models in wind tunnel 
tests [4], This is a demanding application because the distances between the magnets and 
the model are large, but the forces required may be small. 

More recently, with the development of solid state power electronics and advances in 
controls, more attention has been devoted to the possible applications of magnetic 
suspension to industrial and laboratory machinery where large forces may be involved. 
Nikolajsen, et. al. [5] reported on the use of an electromagnetic damping device for 
controlling vibration in a flexible transmission shaft. Schweitzer, et. al. [6] considered the 
application of magnetic bearings to vibration control of pumps and centrifuges, and 
discussed the merits of centralized versus decentralized control [7]. The use of magnetic 
bearings in a flexible rotor system was also considered [8]. 

A number of papers have been published beginning in the early 1980's on various 
aspects of the control of shaft vibration and suspension by magnetic forces. 

Allaire, et. al. performed theoretical studies of the effects of using a feedback actuator 
on the unbalance response of a single mass rotor on rigid supports [9], and on flexible 
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supports [10]. The actuator was placed at the mass location and was represented by 
feedback with gains proportional to shaft displacement and shaft velocity. It was found 
that proportional feedback could be used to alter the critical speeds of the system over wide 
ranges, and that derivative feedback could be used to change the amplitudes of vibration. 
Combinations of proportional and derivative feedback significantly altered the system 
characteristics in terms of both critical speeds and amplitude of response. 

An experimental test apparatus for applying feedback control to a multimass rotor at 
the bearing locations was constructed by Heinzmann [11]. The rig used actuators made 
from the moving voice coils of loudspeakers, and the force was applied directly through 
mechanical links attached to ball bearings on the shaft. Significant effects on the critical 
speeds of the system were achieved by feedback control. 

Kelm [12] computed the linearized stiffness and damping coefficients for a four- 
magnet bearing, and measured the coefficients in an experimental bearing for a two-inch 
diameter shaft. Precise agreement between predicted and measured values was not 
achieved. 

Connor and Tichy [13] have proposed an eddy current bearing that would generate 
repulsive rather that attractive forces by inducing currents in the rotor. 

Chen and Darlow [14] tested a magnetic bearing constructed by modifying an 
induction motor stator and evaluated the effectiveness of two schemes for estimating 
velocity and acceleration in the feedback control loop. 

Walowit, et. al. [15] and Albrecht et. al. [16] analyzed and tested a magnetic thrust 
bearing. Their analysis and experiment involved transverse misalignment of the plane 
surface gaps but no angular misalignment. 

Keith, et al. [17] have examined several aspects of proportional-derivative control 
using a digital controller, and Maslen, et al. [18] consider some of the performance 
limitations of active bearings. 

Papers contained in the proceedings of the first significant international gathering of 
researchers in the field of magnetic suspension of machine elements [19] address 
applications, control, identification of parameters and other aspects of magnetic bearings 
[20-22], as well as applications in space [23-24], 

1.1.3 Developments During this Project 

According to literature from magnetic bearing manufacturer S2M, as of 1991 a total of 
more than 440,000 hours of operation had been accumulated by machines equipped with 
the company's active magnetic bearings [25]. The 96 individual machines span a wide 



7 


range of sizes, from blowers in the 5 to 200 kW range up to industrial compressors in the 
25,000 kW range. 

Along with increased industrial application of magnetic bearings came significant 
new research. The proceedings of the Second International Symposium on Magnetic 
Bearings [26] contain 53 technical papers on various aspects of magnetic bearings, by 
authors from 12 countries. 

These papers address several areas of application, including momentum wheels for 
energy storage [27, 28], electrospindles for boring, grinding and milling operations [29, 

30, 31], as well as suspension of large industrial machine rotors such as those of boiler 
feed pumps [32], pipeline compressors [33], and nuclear circulating pumps [34], 

There continued to be a growing interest in control aspects, with papers devoted to 
digital control [35, 36, 37] and amplifier design [38, 39]. Several philosophies of control 
were examined, including centralized vs. decentralized control [40], automatic balancing 
[41], and modal control [42, 43]. One method of approaching linearity in magnetic 
suspension systems is to apply large bias currents, upon which are superposed the control 
currents. Higuchi et al. [37], however, used a digital control scheme to effect a 
linearization of the magnetic bearing properties without using large bias fluxes. 

Herzog and Bleuler [44] proposed the use of H°° control to achieve required 
stiffness over wide bandwidth, and Fujita et al. [45], seeking a robust control design, 
implemented H°° control using a commercial digital signal processor. Experiments 
indicated that the system was highly stable when subjected to step disturbances. 

Ueyama and Fujimoto [46] physically measured the iron losses due to hysteresis and 
eddy currents by monitoring the coastdown of a rotor suspended in magnetic bearings, at 
different values of coil current, and propose an empirical equation to represent these losses. 

Zhang et al. [47] describe a magnetic bearing application in which the rotor is a thin 
flexible shell, and discuss the advantages of individual magnet control versus control of 
opposing magnet pairs. They conclude that improved damping is possible using the 
individual magnet control. The authors speculate that the method will also be advantageous 
in suspending travelling metal sheets. 

Stability of a suspended rotor was considered by Chen et al. [48], but as in previous 
such analyses, the representation of the magnetic forces is based on a linearized model. 

Of particular interest in the context of the present project, Satoh et al. [49] examined 
a self-excited vibration of a suspended rotor in a flexible structure. The authors concluded 
that interactions between the mechanical structure and nonlinearity of the electromagnets led 
to a vibration with two frequency components. 
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A recent meeting devoted primarily to magnetic bearings was ROMAG’91, 
organized by the University of Virginia and held in Alexandria, VA in March 1991. In 
addition to considering applications in turbomachinery, some presentations also dealt with 
use of magnetic suspension in vibration isolation, particularly in applications related to 
space experimentation [50, 51], although one paper presented a digitally controlled 
magnetic suspension and vibration isolation system for optical tables [52], 

With regard to magnetic bearings for turbomachinery, a number of applications were 
discussed, ranging from canned pumps to gas turbine engines [53] and rocket engine 
turbopumps [54], 

Again, considerable emphasis was placed on control aspects, with papers devoted to 
the effects of sensor location [55], effects of amplifier design [56] and the general 
controllability of flexible rotors [57], 

Subsequent meetings have explored a number of these aspects in greater detail. These 
include the Third International Symposium on Magnetic Bearings [58], and Mag ’ 93 [59], 
While the papers in these meetings address progressively more sophisticated control 
strategies, in much of the work presented, variations on a one-dimensional force model are 
used. 
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2. TECHNICAL FINDINGS: Analytical/Numerical Modeling 

The objective of the modelling is to calculate the force exerted by the magnets on a 
journal in the case of steady currents through the coils. The techniques needed for this 
computation can also be applied to calculation of force in the dynamic case if the problem is 
assumed quasi-static in a magnetic sense. It is expected that this will be appropriate in 
most magnetic bearing applications, since the principal requirement for this assumption is 
that the frequencies of current and field variations do not approach radio frequencies. 

Some correction may be necessary to account for eddy current effects, which are neglected 
in the present work, if these methods are applied to the rotating shaft case. 

The results of this section are also described in the Ph.D. thesis of Xia [60], 

2.1 Fundamental Principle 

The principle of force calculation is that the force component in a given direction is 
equal to the negative of the rate of energy change with respect to that coordinate, that is, 

F x = -— (2.1.1a) 

8x 

F y = -^ (2.1.1b) 

5y 

where the energy U is the energy associated with magnetic flux density contained in the 
magnetic circuit 

U = U B HdV (2.1.2) 

z Jv 

where B is the magnetization and H is the magnetomotive force. If the linear 
approximation is made that B = p H, then this may be written. 



Development of the force model proceeded in two stages. The first method that 
was developed considered only the energy in and near the air gaps, approximating the 
metals as infinitely permeable. This method is referred to subsequently as the air gap 
method. It does account for nonuniform gap geometry as well as nonuniform distribution 
of flux within the gaps. The second method, referred to below as the full magnet method, 
includes the energy in the metal of the magnet and a portion of the rotor as well as that in 
the gaps and nearby air regions. Nonlinear magnetization functions can be considered as 
long as they are single-valued. 
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Both methods rely on two-dimensional finite element calculations of flux 
distributions. Computations are performed for one magnet at a time, and interactions of 
flux loops of individual magnets are neglected. 

2.2 Air Gap Method 

A computer program was written that calculates the force exerted on the journal by 
a magnet having a steady current in its coils. The force is found by calculating the energy 
stored in the air gaps between the magnet and the journal, then performing a numerical 
perturbation to obtain a central difference of the energy change per unit position change. 
This gives the force in the direction of the perturbation. 

The force for a magnet at an arbitrary location can be calculated. The calculation 
includes the following assumptions: 

i . The permeability of the metal is infinite compared to that of the gaps, which is 
assumed equal to that of free space. This implies that all the energy is stored in 
the gaps. 

ii. There is no flux leakage, but expansion of the flux lines beyond the gap edges 
is allowed. 

iii. The coil current, therefore the MMF, is constant over a perturbation. 

In an isotropic domain not containing currents, where time variations are only of 
low frequency, the magnetic field can be represented as the gradient of a scalar field 4>(x,y). 


The energy contained in the domain is given by Equation (2.3) above where the 
flux density B is given by 


B = - V<|> 

(2.2.1) 

and the potential <J) satisfies the governing equation 


v 2 <t> = o 

(2.2.2) 

with the boundary conditions 


— = 0 on free boundaries 
3n 

(2.2.3) 

and, because of assumption (i) above 


<J) = d>i on pole face 1 


<J) = d >2 on pole face 2 

(2.2.4) 

4> = 0 on journal surface 



as shown in Figure 2.2.1. 
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Initially the boundary values of <J> on the pole faces, Oj and O2 > are not known, 
but must be determined in relation to the datum of <>=0 on the journal surface. The problem 
is made tractable by the fact that the governing equation is a linear one, so that the values of 
<|) internally are determined within a multiplicative constant even for an arbitrary choice of 

boundary condition values. The fact that the flux must be the same through the two gaps 
allows the ratio 

K= ^ (22 - 5) 

Oi 

to be determined. Then the fact that the difference between the two potentials is the 
magnetomotive force, 

°i-^2 = C (2.2.6) 

allows determination of the actual surface potentials. The variables are nondimensionalized 
so that O] = 1. The procedure is as follows: 


a. At an unperturbed position A, start with Oi = 1, 0 2 * = 1. An asterisk 

represents an initial guess or a calculated value based on an initial guess. 
Later the ratio 


K - ^2 _ 02 
01 02 

will be found. 


(2.2.7) 


b . Solve for the distribution of <j) in each gap based on these boundary 
conditions: <j>i, <|>2*. Note that 4*2 = k<{> 2*. 

c. Calculate the resulting flux density distributions and the energy stored in 
each gap, plus the flux through each pole face 




( 2 . 2 . 8 ) 



dA = 1 


K 


Y2 


(2.2.9) 


d. Since the actual dimensionless fluxes are the same magnitude ( yj = —y2 ) 
k is uniquely determined as the ratio 



* 


Y2 


( 2 . 2 . 10 ) 
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e. The potential distribution in gap 2 is given by <J >2 = K(|) 2 * and the flux 
density is related to the * distribution by the same ratio. The energy is 
therefore given by 

G2 = k i. 2 C2 ( 2 . 2 . 11 ) 

f . Note that the imposed mmf is the difference between the pole potentials 

£ = <!>, -<D 2 = (1 -k)Oi (2.2.12) 


or, since 0] = 1, 


?= a-K) 


(2.2.13) 


(If dimensional values are needed, the factor ( 1 - k ) is also the ratio 
between mmf in amp-turns and the dimensional potential on pole face 1). 
For now, use the nondimensional £. 

g. At a perturbed position B, start with Oi* = 1, 02** = 1. The potential 
on pole 1 is now a guess, since only the mmf was maintained constant 
during the perturbation. 02 ** is a "double" guess because it is based on 
01 * 

h . Use a procedure analogous to steps b through f to find the ratio 


* 

01 


(2.2.14) 


i . Use the calculated mmf £ from step f above, which is still equal to the 
difference between the pole potentials, to write 


and define 





(2.2.15) 


(2.2.16) 


j . Using the logic of step e 


G] = a z Gj 


and 


(2.2.17) 



14 


a 2 = oc 2 X 2 G2* (2.2.18) 

k . Stored energies have now been calculated at position A and position B. A 
forward difference analogue to the force in the direction from A to B is 
therefore 

Fab = ( gl+g2 ^~( a ' + ° 2 ) A (2.2.19) 

Ax AB 

l. To return to dimensional values, use the factor (1 -k) from step f above. 

Although the description above uses a forward difference, a better result is obtained 
using a central difference, which is the method actually employed. This requires 4 
perturbations to find the vector components of force in the x and y directions. 

Although the calculation of forces with the inclusion of three dimensional effects, 
flux leakage and hysteresis will involve significantly more computations, the overall 
approach should be the same as that used above. It may be necessary to use a vector 
magnetic potential, and the assignment of boundary conditions will be considerably more 
complex, however. 

2.3 Air Gap Method: Computer Program 

The algorithm above is embodied in a FORTRAN computer program, GAPFOR1, 
which uses the finite element method for calculating the magnetic potential in two 
dimensions. For a given journal position the program calculates the gap height as a 
function of angular location and generates a finite element mesh for each gap. Flux 
fringing is allowed by extending the finite element domain beyond the edges of each pole 
face. Then the journal position is perturbed four times, first with positive dx and negative 
dx, then with positive dy and negative dy. At each step the mesh is regenerated and the 
energies are recalculated. 

To achieve rapid computational speed and efficiency, a dedicated finite element 
program was written for this application. It includes a grid generation routine as well as a 
banded gauss elimination solver for the assembled equations. A listing of the algorithm is 
given in Appendix A. 
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2.4 Air Gap Method: Results 

2.4.1 Test Case 1: parallel surfaces 

The algorithm was tested by calculating the force in the case of a magnet and part 
having parallel faces, with no flux fringing allowed, shown in Figure 2.4.1. In this case 
an analytical expression approximates the force per unit depth of pole face as 


F = 


_ UgAgNV 


(2.4.1) 


where F = force 

(ig = magnetic permeability of the gap 
Ag = area of the gap 

N = number of magnet coils 

i = current 

h = gap height 

The computer program was run for the slightly different case of an annular clearance 
between a shaft and a magnet, corresponding to the case of a centered journal. For the 
sample case the equation gives F = 22.4 N, while the program predicts F = 20.9 N. 


2.4.2 Test Case 2: Effects of element size and fringing 

The algorithm was used to calculate the force from a single magnet acting on a 
journal, as in the experimental apparatus. The dimensions are given in Table 2.4.1. The 
effects of variations in element size were examined along with the effects of allowing 
fringing to occur by extending the domain of solution circumferentially beyond the ends of 
the pole faces as shown in Figure 2.4.2. Table 2.4.2 shows the results of these variations. 
The column A displays results without fringing, while column B shows results with 
fringing allowed in a domain extended 10% of the width of the pole face to either side. 

The results indicate that without fringing, the effect of decreasing the element size is small, 
but when fringing is to be accounted for, the element size is a significant parameter. The 
results of column B suggest that when fringing is allowed, the predicted force is smaller 
than when fringing is not considered. This might be expected, since fringing decreases the 
average flux density by increasing the volume of the energy storage area. Since the energy 
is related to the square of the flux density, an overall decrease in stored energy and in force 
seems appropriate. 


2.4.3 Forces from one magnet of a bearing 

The computer program has been used to predict the forces from one magnet acting 
on the journal at various positions of the journal within the clearance space. Half of the 
entire clearance space is mapped, since all positions of the journal with respect to a single 
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current: 

turn: 

magnetomotive 


i=l A 

N=200+200 

force: mmf=400 turn-A 


Analytical solution for a flat surface 
by approximate equation 
Fy=22.4 N (5.04 lb) 


Numerical solution for annular clearance 
Fy=20.9 N (4.70 lb) 


Figure 2.4. 1 Comparison of numerical and approximate analytical solutions for 
& magnet. 


one 
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GRID 


n=10 circumferential divisions 
m=5 radial divisions 

AX=0.2X 

( to allow fringing) 


rotor 


Figure 2.4.2 Sample grid (radial dimension exaggerated). 


Pole depth 
Pole width 
Gap height 
Anglebetween poles 
MMF 


10.1 mm 0.750 in 
13.6 mm 0.534 in 
0.76 mm 0.03 in 
40° 

400 A-T 


Table 2.4. 1 Parameters for sample calculations. 
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magnet can be represented in terms of positions in this half space. Figures 2.4.3 and 2.4.4 
show maps of force versus x,y position. The magnet is the upper vertical magnet, and a 
steady 

Attraction Force from One magnet Using Finite Simplex Method 

Design Parameters: 

R = 38.1 mm 
c = 0.76 mm 
L = 19.05 mm 


01 =60°, 02 = 80°- a = 20° 

A B 


Case 

m 

n 

elements 

global matrix 

No Fringing 
Fy 

With Fringing 
Fy 

1 

8 

20 

320 

189 x 11 

20.925 

21.220 

2 

8 

70 

1120 

639 x 11 

20.927 

20.635 

3 

8 

150 

2400 

1359 x 11 

20.927 

19.881 


Table 2.4.2 Effects of element size and fringing 


current of 1 ampere through the coils is used. The dimensions and other parameters are the 
same as those of the experimental apparatus described below. The figure indicates that the 
force in the y-direction varies between 0 and 132 N as the journal is moved along the y 
axis. When the journal is also given an x-direction eccentricity, the y-force decreases 
significantly. Except at x=0, there is also a small x component to the force, shown in 
Figure 2.4.4. 

In a subsequent section the predicted forces are compared with those measured in an 
experimental apparatus. 

2.5 Full Magnet Method 

Unlike the previous method the present section considers the magnetic flux within 
the metals in addition to that in the air gaps. This allows the examination of effects such as 
local magnetic saturation of the materials and residual magnetization. This approach 
presents two categories of difficult problems, however. The first category arises from 
consideration of finite permeability, which in the general case is a nonlinear and 
multivalued function of field intensity. The second is related to boundary conditions on 
magnetic field quantities, and a third concerns the source, or current density, term of the 



19 


SUSPENSION FORCE FY(N) 

b r7i* INCH C = 0.03 INCH UUF-2004200 T-A 



Figure 2.4.3 


Vertical force of attraction from upper magnet with 1 A current, by numerical 
calculation. 
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SUSPENSION FORCE FX(N) 

R = 1 5 INCH C = 0.03 INCH HMF=200+200 T-A 
XB=X/C YB=Y/C 



Figure 2.4.4 Horizontal force of attraction from upper magnei with 1 A current, by 
numerical calculation. 


governing equations. First the forms of the differential equations will be presented and 
then approximations and assumptions will be introduced to simplify the equations. 
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2.5.1 Differential Equations in 3-D 

In the general three-dimensional case, the flux density B, the field H and the 
current density J are related by 

V x H = J (2.5.1) 

B = |iH (2.5.2) 

VB*=0 (2.5.3) 

where 

(X = Ififtt) (2.5.4) 

Assuming that a solution for B in all parts of the domain can be found, the method 
of force calculation described in Section 2. 1 can be applied. 

2.5.2 Two-dimensional Equations 

For the next phase of the analysis, the solution domain is simplified from three 
dimensions to two dimensions. Successful finite element solutions for magnetic flux have 
been obtained in two dimensions (Chari and Silvester [61]) for cases of single valued 
permeability by making use of Equation (2.5.3) to write the flux density, or magnetization, 
vector as the curl of a vector magnetic potential 

B = V x A (2.5.5) 


This equation is valid in three dimensions, but is more easily applied if the magnets and 
rotor are treated as infinitely long and the current is assumed to pass only in the coils of the 
magnets. Under these approximations both the vector potential and the current density 
have only one component (z), and Equations (2.5.1), (2.5.2) and (2.5.5) can be combined 
to write 


a 2 A 


a 2 A 


= |XJ 


(2.5.6) 


3x 2 3y 2 

where A and J are magnitudes of the corresponding vector quantities. For given |l and J 
distributions and appropriate boundary conditions, this Poisson's equation can be solved 
by finite difference or finite element methods. In the present work the current density is 
assumed uniform within the coil windings and zero elsewhere. The coils are treated as 
isotropic solids, as shown in Figure 2.5.1. 




Current density J 
Current density -J 


Figure 2.5. 1 Modelling of coils with uniform current density 
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At this point the two most difficult problems arise; determination of permeability 
and assignment of boundary conditions. 

2.5.3 Permeability 

In general, the permeability of a ferromagnetic material is a nonlinear, multivalued 
function of field intensity and, through the history of the field strength, of time. The 
magnetization of the material is often represented by the hysteresis loop (a) of Figure 
2.5.2, adapted from Cullity [62] which shows the relationship between H and B for a 
particular time variation of H, namely a cyclic completely reversed variation that is 
sufficiently strong to cause the material to be saturated alternately in both directions. 
Although this figure gives some qualitative insight into the material's behavior, it does not 
fully characterize the response of a magnet to other types of time varying excitations. In 
fact, for an excitation H that does not fully saturate the material, the curve traced by the B 
function might take one of several other forms shown in Figure 2.5.2, depending on the 
material and the range of H. 

For analytical purposes, it is most convenient to assume a linear variation of B with 
H, or a constant permeability ( (a) of Figure 2.5.3). For this assumption the solution of 
Equation (2.5.6) is straightforward and obtainable by a direct method. Next in complexity 
is the consideration of B as a nonlinear but single valued function of H, as in (b). An 
iterative method is now required for the solution. In addition, the calculation of the energy 
in the magnetic field. Equation (2.1.2), requires integration using the actual magnetization 
function. The most general case, that where B can take on an infinity of values, 
depending on the history of H, is not considered in the present work. Therefore, in this 
report calculations are limited to single-valued functions of B vs. H. 

2.5.4 Boundary Conditions 

Far away from the magnets and rotor it is reasonable to assume that the magnetic 
flux intensity B is zero, which implies that 

3A 3A 

ST "*" 0 <2 ' 5 ' 7) 

It is feasible to extend the solution domain far enough to approximate this condition. 
Numerical studies of the effects of domain size were made as part of the analysis, and it 
was noted that the penetration of magnetic flux into the rotor is limited. A sample 
discretization is presented in Figure 2.5.4, where the boundary conditions given by 
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Figure 2.5.2 Possible B-H loops in a real ferromagnetic material. 
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Figure 2.5.4. Example of modified direct discretization of portion of domain. 


equation (2.5.7) are applied af , • 

2 -5-5 Discretization 
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I - 1 A 
x * 0.0 in. 
y - -0.024 in. 

Extra boundary 2 0 % 

10 % with the fine elements 
within the pole area 


m x n Force (N) 


n = 20 


2X5 

7.334 

2 X 15 

6.557 

2 X 20 

6.361 

2 X 30 

5-961 

4X5 

7.397 

4X 15 

6.654 

4 X 20 

6.729 

4X 30 

6.514 

6X5 

7.409 

6 X 15 

6.933 

6X 20 

6.645 

6X 30 

6.696 

6X5 

7.413 

6X15 

6.964 

6 X 20 

6.694 

6X 30 

6.763 

2X5 

7.334 

4X5 

7.397 

6X5 

7.409 

6X5 

7.413 

2 X 30 

5-961 

4X 30 

6.514 

6X 30 

6.696 

6X 30 

6.763 



Difference between the largest 
and smallest forces 


1-353 


0.663 


0.713 


0.63 


0.079 


0.602 


Table 2.5.1 Divisions and parametric study of grid size effects in gap region. 


Figure 2.5.5. Example of variable density grid after Cavendish. 
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position followed by calculation of the forces based on changes in energy using a central 
difference approach, including the energy change due to flux distribution changes during 
the perturbation. In the case of larger domains having nonlinear permeability the 
computation time for this method would be large. The second method, which was actually 
adopted, holds the flux distribution constant and calculates the energy change due to the 
area changes of all the elements that are distorted during a perturbation. This method is 
much faster than the full numerical perturbation scheme, and has been shown [67] to have 
high accuracy. 

2.6 Results of Linear Calculations 

In Section 2.2, a method was described to calculate forces using a linear method, in 
which the air gaps only are treated and the flux distribution is calculated by the Laplace 
equation for the scalar magnetic potential. Results of calculations using this method were 
presented in Section 2.4. In Section 2.5, the linear method was extended to include 
regions of differing permeabilities, using the Poisson equation for flux distribution. The 
present section presents results of these calculations, along with experimental 
measurements. For a description of the experimental apparatus and methods, refer to 
Section 3. In that section, some of the calculations presented here will be shown again. 
The work presented in this section is also described in the paper "Determination of Forces 
in a Magnetic Bearing Actuator: Numerical computation with Comparison to Experiment," 
by Knight, Xia, McCaul and Hacker [68]. Only the Conclusion section of the paper is 
reiterated here. 

Conclusion (of Reference [68]) 

Calculated and measured forces in a magnetic journal bearing actuator 
are presented. The calculations are based on two-dimensional finite element 
solutions of the magnetic flux distribution in both metals and free space. The 
measurements were made in an apparatus designed for direct force 
measurement by strain gage transducer assemblies supporting a non-rotating 
journal. 

Comparison of numerical calculations with one-dimensional magnetic 
circuit theory indicates that as the gaps are made non-uniform by the approach 
of the journal to the magnet, two dimensional effects become significant and 
the two methods predict different forces. At relative permeabilities above 
1()4, changes in permeability of the metal have little effect, but at lower 
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permeabilities the available force decreases dramatically with decreasing 
permeability. 

Also predicted is that the effect of finite metal permeability is more 
strongly felt at small gaps than at large gaps. 

The calculated principal attractive forces agree well with the measured 
forces when a relative permeability p r = 500 is used, corresponding to highly 
saturated material. The measured normal forces, however, are higher than the 
calculated values even when a high permeability is used. 

It seems reasonable that the permeability distribution in the metal is 
non-uniform. Future work is planned in which distributions of permeability 
will be examined. 

After the submission of this paper, the numerical method was extended to model 
nonlinear distributions of permeability. 

2.7 Nonlinear Force Calculation 

An algorithm was developed to calculate the force exerted on a rotor by a magnet, 
considering the effects of a nonlinear magnetization characteristic for the rotor and magnet 
material. It uses the finite element method to solve the equation for vector magnetic 
potential in two dimensions. The force calculation part of the algorithm is based on the fast 
solution method proposed by Coulomb [67], There are three primary operations involved 
in the force calculation: (a) modelling of the magnetization curve of the magnet and rotor 
material, (b) iteration for the distribution of vector magnetic potential consistent with the 
nonlinear permeability, and (c) application of the force calculation algorithm. These 
operations are outlined briefly below, but more complete descriptions of the methods are 
given in Appendix C. 

This work is also described in a paper, "Forces in Magnetic Journal Bearings: 
Nonlinear Computation and Experimental Measurement," by Knight, Xia, and McCaul 
[69], presented at the Third International Symposium on Magnetic Bearings, Alexandria, 
VA, July 1992, and contained in the proceedings of that meeting. 

2.7. 1 Modelling of Magnetization Curve 

For most of the calculations presented here, the magnetization function for silicon steel 
[62] has been used. Some calculations were also performed using an arbitrarily chosen 
function having sharp discontinuity in slope, to assess the effects of abrupt saturation. 



32 


The magnetization function of the steel is nonlinear but single-valued; that is it does 
not exhibit hysteresis. The function is represented by tabular data and is approximated by a 
cubic spline interpolation. At field intensities higher than 1200 A t/m, the slope of the 
magnetization function is assumed to be the permeability of free space, p<). Figure 2.7. 1 
shows the actual magnetization data and the approximation. 

For numerical calculations, a more useful representation of the magnetization function 
is that of Figure 2.7.2. The reluctivity of the material, H/B, or 1/p, is plotted versus the 
square of the flux density. When this representation is used it is not necessary to calculate 
the field intensity at each location for every iteration, but only the flux distribution. 


2.7.2 Calculation of Flux Distribution 


The distribution of scalar magnetic potential, leading to the distribution of flux density, 
is calculated by the finite element method. The equation that models the potential is the 
nonlinear Poissons equation 


/±§ a 


a 1 1 dA 


ariirarr^iir^h 1 (Z7J> 

where A is the magnitude of the vector magnetic potential, which in the two-dimensional 
case has only one component , normal to the plane of the solution region. The flux density 
is related to the potential by 


B = V x A. 


(2.7.2) 


This relationship allows a convenient representation of flux density, since it implies that 
contours of constant A are also lines parallel to B. 


The source term J, current density, appears in those elements comprising the cross 
sections of the coils. The value of the total ampere-turns is divided by the nominal cross- 
section to arrive at this current density. 

An iterative method is used to obtain a distribution that is consistent with the nonlinear 
magnetization function. The procedure is that recommended by Silvester [70], in which a 
Newton-Raphson iteration is applied to determination of the reluctivity. An initial 
approximation to the potential is made, then updated based on successive solutions of the 
Poissons equation for incremental changes in the A field that result from refinement of the 
reluctivity distribution. 


When the flux distribution has been determined, the calculation of forces is performed 
using the method of Coulomb [67], in which only the energy changes in the distorted 
elements are considered during a virtual displacement. The method allows the force to be 
determined without multiple solutions for the flux distribution 
Appendix C describes the numerical methods in more detail. 
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2.7.3 Results of Calculations 

Calculations were performed based on the geometry of Figure 2.7.3, corresponding to 
the first experimental apparatus and the measurements described in [71], The magnet 
under consideration is an upper vertical magnet, so forces in the y-direction are the 
principal forces, and forces in the x-direction are the normal forces. Also plotted in the 
figures are the results of the linear calculation described in previous reports and in [71], 

The effect of saturation on the force is seen in Figure 2.7.4, which shows the attractive 
principal force as a function of the coil current, when the shaft is in a centered position with 
respect to the magnet pole faces. The gap between shaft and magnet poles is therefore 
constant at 0.03 inch. The dimensionless force is seen to increase with current, and below 
a current of 2.5 A (corresponding to 1000 A t) the result of the nonlinear calculation is the 
same as that of the linear calculation. Above this value the force continues to increase, but 
at a much smaller rate than predicted by linear theory. 

At current levels higher than 3 A the magnet material experiences saturation near the 
inner corners of the intersection between the pole legs and the magnet outer arc. As the 
current level is increased, the area of saturation expands across the cross section of the 
legs. Figure 2.7.5 shows those elements that have been saturated for the case of i = 3.5 A. 
At this level of MMF the area of saturation encompasses a complete layer of elements 
spanning the cross-section. For purposes of this plot, saturation is defined to correspond 
to a flux density of 1.4 T. At this point the slope of the magnetization function is assumed 
to be that of free space, so above this level of flux density the force can continue to increase 
with current, as indicated by Figure 2.7.4, but at a much slower rate. 

For a given MMF the magnet may also experience saturation when the shaft is moved 
closer to the magnet. Such a displacement decreases the overall reluctance of the magnetic 
circuit by closing the gaps, and changes the gap shape as well. Figures 2.7.6 to 2.7.8 
show the increase in number of saturated elements when the shaft is moved toward the 
magnet, for the constant current i = 2.0 A. Figure 2.7.6 corresponds to a shaft eccentricity 
of (X,Y) = (0, 0.5), which denotes a position on the magnet's axis of symmetry, half the 
distance from the center to the maximum possible eccentricity. There are two areas where 
elements are saturated; the inner corners of the horseshoe, and the part of the shaft near the 
inner edges of the pole faces. These edges are the points of closest proximity between the 
poles and the rotor. As the shaft is moved closer to the magnet the areas of saturation 
enlarge. At an eccentricity of 0.6, Figure 2.7.7, the upper ends of the pole legs have been 
completely saturated, and the area of saturation at the rotor surface has expanded. As the 
eccentricity is further increased to 0.7, Figure 2.7.8 shows the saturation areas continuing 
to expand. The contour plot of Figure 2.7.9 reflects the saturation pattern. Comparison of 
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Figure 2.7.9 and Figure 2.7.10, which is the potential distribution obtained by a linear 
solution, shows how the flux distribution has changed in order for the flux lines to 
maintain a minimum curvature and to follow the easiest path, while minimizing local 
concentrations. The effect on the force is illustrated in Figure 2.7. 1 1 , where the nonlinear 
calculation is compared with the linear solution using a relative permeability of 5570 
(corresponding to that of silicon steel at very low field intensity). Above an eccentricity of 
0.4, the force continues to increase, but at a much lower rate than predicted by linear 
theory. 

Asymmetry in the distribution of saturation develops when the rotor is given an 
eccentricity away from the magnet's symmetry axis. Figure 2.7.12 shows the saturation 
pattern when the shaft is moved to the right, to a position (0.45,0.7), a large eccentricity. 
The rotor near the inner corner of the left leg is saturated, as well as almost an entire layer 
of elements near the right leg. The saturation region at the upper ends of the legs has also 
changed slightly from that of Figure 2.7.8. Figure 2.7.13a shows the potential distribution 
for this case. Comparison with Figure 2.7.13b, which is the potential distribution obtained 
by a linear solution, illustrates the effect of saturation in excluding some of the flux from 
the corners and increasing the fringing at the poles. 

The force in the normal direction is plotted in Figure 2.7. 14 for one value of off-axis 
eccentricity, as a function of the y-position (along the symmetry axis). This force also is 
predicted to deviate from the linear theory above a y-displacement of 50 % of the clearance. 

The ratio of normal force to principal force for this same normal eccentricity is plotted 
in Figure 2.7. 1 5. Curves are shown for the nonlinear theory as well as for the linear 
theory at two different relative permeabilities, in addition to the experimental results. 
Although the nonlinear theory does not reflect the magnitudes of the ratio very accurately, 
the trend is appropriate: the nonlinear theory does indicate a very slight increase in this 
ratio as the eccentricity is made larger. The magnitude of the difference, however, is so 
small that it may not be significant in view of the numerical solution method. 





Principal force (nondimensional) 


x/c=0.0 



y/c - principal coordinate 


Figure 2.7.11 Principal force as a function of the principal coordinate. 
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Figure 2.7.13b Linear magnetic potential contours at large normal eccentricity. 



Normal fores (nondimensional) 


x = 0.45 



= 2A 


1 .0 


y/c - principal coordinate 

Figure 2.7. 14 Normal force at normal eccentricity of 0.45 . 





Force ratio 
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x/c - 0.45 



Figure 2.7. 15 Ratio of normal to principal force at large normal eccentricity. 
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2.8 Effects of Uncertainties and Property Variations 

In earlier results, numerical calculations of forces using the nominal geometry of 
the experimental apparatus did not predict accurately the magnitudes of the normal forces. 
The measurements were in all cases considerably larger than predicted by calculation. 
Recent calculations have attempted to address the issues of uncertainty in the pole face 
geometry on the forces. The nominal geometry of the apparatus is shown in Figure 2.8.1, 
along with one possible type of geometric error. Suppose that each of the pole faces, while 
still consisting of a circular arc, is rotated a small amount from its nominal orientation. 

This rotation is distinct from the angular uncertainty described in section 2.2 above. This 
type of error, or an error of similar magnitude, might result from tolerances in machining, 
but would be unlikely to result from assembly errors. 

Figures 2.8.2, 2.8.3 and 2.8.4 show the results of several series of calculations. 
Each figure is based on the same set of shaft positions and presents the force from one 
magnet when the shaft has a large normal eccentricity, as a function of position on the 
principal coordinate. These data are most easily compared with data from the first force 
measurement apparatus, because of the sequence of measurements. The four curves 
illustrate the effects of two different magnitudes of error in pole face orientation, 0.5° and 
1.0°. These correspond to movement of the outer corner of each pole face a distance of 
0.004 inch or 0.009 inch toward the shaft. A change of 0.5° therefore represents 15 % of 
the radial clearance. The effects of these changes on the calculated forces is shown for the 
case using the nominal magnetization function for the material, with a saturation flux 
density of 1 .4 T (curves 1 and 4), and for the case using a saturation flux density reduced 

by 20 %, to 1 . 14 T. Asa reference, the result for the linear calculation with no error in 
pole face geometry is included. 

It is seen that the effect of this geometry change is to increase both the principal and 
normal forces, with a larger angular change causing larger forces as long as the nominal 
saturation flux density is maintained. The ratio of normal to principal force also increases. 
For an angular error of 1°, the ratio of forces has approached the ratio that was measured 
When the calculation is performed using the lower value of saturation flux density, 
however, the result is similar up to the point where saturation is felt, then the ratio between 
the forces decreases with increasing principal coordinate. It is reasonable that a 

combination of geometric error and saturation flux density level can produce the force ratio 
observed in the experiment. 

experimentally. In contrast to the experiment, however, the ratio shown in Figure 2.8.4 
increases with principal coordinate, while the measurement does not indicate this trend. 



x/c (principal coordinate) 


Figure 2.8.2 Effect of rotation of pole faces on calculated principal force. 





x/c (principal coordinate) 


Figure 2.8.3 Effect of rotation of pole faces on calculated normal force. 
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2.9 Conclusions (Analytical/Numerical) 

Principal conclusions drawn from this phase of the work are- tn tha, ih, a , „ 

the force tn the norma! direction, (2) that the normal foree measured experimentally is 
evera times as large as the magnitude predicted at present by either linear or nonlinear 

to Z' ,'f ' he lrc " d ° f n ° nlinear,he0 ^ P redi « 'arger nonnal forces in relation 
to principal forces is appropriate. 
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3 TECHNICAL FINDINGS: Experiments 

Two different types of experimental apparatus were constructed to measure drrectly 

the forces exerted between magnets and shaft in the nontotatmg case. 

The first apparatus was made using a solid dtsk to approximate the rotor, and so 
core magnets. Zs apparatus made use of stratn gage instrumented support arms for the 

disk in order to measure the force directly. 

The second apparatus was made with laminated disk and magnets, and used 

orinciple of a calibrated deflecting beam to measure the forces indirectly. 

^ Measurements of force on a stationary, non-rotating shaft were made as functions of 
position and current and the forces have been compared with corresponding numerical 
predictions. Where appropriate, measurements from the two rigs were also compared, 

were found to be consistent with each other. , r701 

These apparatus and results are also described in the M.S. thesis of E. McCau. (72,. 

3 '' Tablel resign parameters of the apparatus for direct force measurement. 
Figures 3.1.1 and 3.1.2 show schematics, and Figure 3.1.3 shows an exploded view . 


Rotor o.d. 

Support shaft o.d. 
Support shaft length 


0.076 m 
0.016 m 
0.15 m 


3.00 in 
0.625 in 
6.0 in 


Magnet i.d. 

Magnet depth 
Shaft clearance (diam) 
Pole width 
Leg length 


0.0777 m 
0.019 m 
1.52 mm 
0.013 m 
0.019 m 


3.06 in 
0.5 in 
0.06 in 
0.5 in 
0.75 in 


Coils 200 turns/leg, #22 copper 

Pole separation angle 40° 

Magnet centerlines 0° 90° 180° 270 

Table 3.1.1 Design parameters of experimental apparatus I 




Figure 3.1.1 Schematic of experimental apparatus I. 
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Figure 3.1.2 Dimensions (inches) of experimental apparatus I. 


0.18 m 
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Figure 3.1.3 Assembly drawing of experimental apparatus I. 
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The magnets and the journal are cut from disks of solid 1020 steel. They are placed 
between two side panels that are laminated from 1/16" aluminum pieces. The magnets are 
held in place by locating pins that are press fitted into through holes. The hole positions 
were located prior to the cutting of the magnets from the solid disk. In this way careful 
control of the radial clearances and angular positions of the magnets was maintained. 

Each magnet is independent and is wound with 400 turns capable of carrying current of 
2.0 A in the steady state. The original design based on a sandwich construction was 
intended to allow flexibility in mounting magnets of different materials. In the steady force 
measurement mode, the rotor is held stationary by pressure from six micrometer heads 
(three on each end) that are in turn held by cantilever arms instrumented with strain guage 
bridges. Thus all mechanical force on the rotor passes through the strain guage arm 
transducers, and ideally all of the force on each micrometer tip is purely radial. In fact, it is 
likely that the transducer arms exert some force in the tangential direction because of static 
friction between the pusher tip and the support disk. Such force would not be sensed by 
the transducers, which are designed to measure only forces that cause bending moments. 
Attempts were made to minimize any frictional force by adding Teflon ball sockets with 
steel spheres between the pusher micrometers and the support disks. 

Despite the difficulties encountered, a number of successful force measurements were 
made and useful conclusions have been drawn. Measurements from the second apparatus 
tended to confirm those of the first. 

3.1.1 Method of Measurement 

Measurement of the force at a given location x,y within the clearance space requires 
several steps: 

i. Establish a datum position relative to the magnet from which the force is to be 
measured. This requires placement of the rotor in contact with the magnet along 
the inner corners of the pole faces. Visual alignment of the rotor is followed by 
application of a small steady current to the magnet to assure contact. The datum 
readings are then taken from the eddy current probes located at 45° to the vertical. 
Ideally, one datum would be sufficient for all positions and all magnets, but in 
reality because of machining tolerances and assembly allowances, the magnet 
pole faces are not located on a perfect circle. Measurements can only be made 
relative to one magnet at a time, therefore, and a separate datum is required for 
each magnet. This datum can be used for all positions relative to this magnet. 
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ii. Place the rotor in the desired position by adjusting the micrometer pushers, 
repeatedly computing the position from the probe readings and correcting as 
necessary. When the rotor is in position, all strain guage arms should be under a 
slight preload. 

iii. Take readings of the strain guage voltages without current through the magnets. 
These will serve as datum values that contain all preloads including the rotor 
weight. 

iv. Apply the desired current to the magnet. Take readings of the strain guage 
voltages and the position probe readings. 

v. Compute the position from the probe readings and the forces from the strain 
guage voltages after subtracting the datum values. 

vi. Apply an alternating current to the magnet coils to remove residual magnetization 
and return to step ii. 

The intent was to automate the entire process of data taking and force calculation by 
using a microcomputer and digital data acquisition. Difficulties with the commercial A/D 
hardware, however, forced the use of manual data taking for this phase of the work. The 
raw outputs from the strain guages and position probes are processed using the same type 
of software as that intended for the automated process, but the data were entered manually. 

3.2 Results of Measurements, Apparatus 1 

Forces were measured at several locations and for several values of steady current. 

The figures referred to below display dimensional data as measured, with forces in 
Newtons plotted against y/c, the eccentricity ratio in the vertical direction. All of the forces 
measured in this apparatus are from the lower vertical magnet, so the vertical forces are in 
the negative y-direction. The eccentricities in the x-direction are all positive. Three 
traverses of the y-direction were made, at x/c positions of approximately 0.0, 0.24, and 
0.45. Assessments of the errors in measurement are not complete; however, it is expected 
that the error in position measurement is no greater than plus or minus 0.05 in y/c and x/c, 
and that the error in force measurement is no greater that plus or minus 5 N. Errors in 
current level control are within 0.1 A. A larger series of measurements that were made 
before the addition of the ball/socket contacts was eventually discarded because the 
measurement error due to friction appeared to be significant. 

The data support some of the anticipated relationships among the position, current and 
force variables but appear to disagree with other aspects of the present theory. Figure 
3.2. 1 shows the vertical force as a function of y/c for several values of current. The force 
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tends to increase roughly as the inverse square of the gap. The magnitudes of the forces, 
however, are considerably lower than those predicted either by the linear finite element 
theory or by the traditional theory based on assumption of uniform gaps, and the ratio 
between measured and predicted forces is not constant. Figure 3.2.2 is a comparison of 
the measured forces with those predicted by the finite element calculation. The results 
indicate that at large gap and/or small current the ratio between the measured and predicted 
forces is about 1.5, but at smaller gaps and/or higher currents this ratio increases, 
eventually exceeding 2.0 for all the three values of current that are plotted. 

Several mechanisms may be operating to cause these discrepancies, including flux 
leakage, non-uniform permeability of the materials and magnetic saturation. Some part of 
the disagreement is likely the result of measurement errors, but the differences appear to be 
significant even after allowing for reasonable experimental error. These disagreements 
reinforce the need for additional work on force calculation. 

The linear finite element theory predicts the existence of forces from a magnet that are 
normal to its axis of symmetry when the rotor is displaced from this symmetry axis, but the 
forces that are measured are considerably stronger than those predicted by calculation. 
Figure 3.2.3 shows the x component of force when the rotor is placed as closely as 
possible on the y-axis. The normal force appears to be somewhat stronger at higher 
current levels but all these forces are small, on the order of 5 % or less of the principal 
force, so it is difficult to attribute much significance to this ratio in view of the experimental 
uncertainty. At higher values of x/c, however, the normal force becomes much more 
significant. Figures 3.2.4 through 3.2.7 show the vertical and horizontal components of 
force when the x/c value is 0.24 or 0.45, and Figure 3.2.8 shows the value of the x force 
as a function of position for several values of x/c while the current is held constant at 1 .0 
A. In general it appears that the normal force increases significantly with increasing x/c, 
and at x/c = 0.24 and 0.45 the horizontal force is about 10 % of the principal force. 

Theory predicts a ratio of about 3 % to 5 %. 
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Figure 3.2. 1 . Vertical force from lower vertical magnet at different values of current. 
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Figure 3.2.2. Measured forces and forces predicted by linear FEM calculation. 
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Figure 3.2.6. Vertical force from lower vertical magnet when x/c = 0.45 
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An aspect of the measured forces that was not anticipated is the lack of degradation of 
principal force as the rotor is moved off the principal axis. Numerical calculations predict a 
significant decrease in the principal force under these conditions, but the measurements do 
not support this prediction. Figure 3.2.9 indicates that within measurement uncertainty 
there are not significant differences in the y-components of force at the three different 
values of x/c. A possible cause is that a self-correcting redistribution of flux along the pole 
faces occurs that allows the force to be maintained. The present theory assumes that the 
magnetic potential along the entire surface of each pole face is uniform. The mechanism of 
potential and flux redistribution should be studied further. 

In summary , the general trends of the measured y-forces agree with the predictions of 
the theory while the magnitudes of forces are somewhat smaller than those predicted. 

Other aspects of theory are not confirmed by the measurements. The measured forces in 
the x direction appear to be significantly larger than those predicted by theory when the 
rotor has an x eccentricity. Also, the y forces do not appear to decrease significantly when 
the rotor is given an eccentricity in the x direction. These effects appear to be significant 
even after considering experimental uncertainty, and both of these phenomena were judged 
to warrant further study. 
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3.3 Measurement Apparatus II 

Measurements of forces exerted by a magnet on a stationary, non-rotating rotor were 
made using a second apparatus that relied on a different measurement principle. 

The magnet cores and the shaft are constructed of 0.014 inch laminations of silicon 
steel M 15. Each magnet is wound with a total of 400 turns of #22 wire, arranged in two 
coils, one on each pole leg. The nominal dimensions of the magnets are the same as those 
of the first apparatus, which was of solid material, but the effective cross section is smaller 
in the new apparatus because of the laminated construction. At present its value has not 
been determined. 

3.3.1 Deflecting Beam Apparatus 

Figure 3.3. 1 is a schematic of the apparatus for force measurement. The shaft is 
clamped at each end in two large pedestals that are fixed to a solid base. The magnets are 
assembled in a retaining shell and the entire magnet assembly is mounted on a slide 
mechanism allowing movement in the horizontal direction. The slide is mounted in turn on 
a laboratory jack that allows the assembly to be moved vertically. Thus the bearing 
assembly can be moved in two directions and positioned accurately with respect to the 
fixed shaft. The relative position of the bearing is measured by four proximity probes 
oriented at 45° to the vertical. These probes are connected to the bearing housing so they 
always measure the relative displacement of the rotor from the center of the bearing 
regardless of the deflection of the rotor support beam. 

When one or more of the magnets is activated, the force causes a deflection of the 
beam from its static position. The components of this deflection in the vertical and 
horizontal directions are measured by a separate set of proximity probes that are connected 
directly to the base of the apparatus. The intention was to place the support beam in the 
pedestals to approximate the perfect clamped-clamped case, so the stiffness of the beam 
would be equal in all directions and could be calculated from simple beam theory. After 
assembly it was found, however, that manufacturing tolerances resulted in unequal 
stiffnesses, so the force vs. deflection relationship was directly calibrated independently in 
both directions. Although the deflections were different in the two directions, the 
relationships were linear over the range required for measurements, so the calibrations 
yield a constant horizontal and a constant vertical stiffness. 




Figure 3.3.1 Schematic of experimental apparatus II. 
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3.3.2 Measurement Method 

Before conducting any force measurements, the location of the bearing center is 
determined by noting the readings of the position probes when the shaft is placed against 
the pole faces of the magnets and interpolating to find the center. Also, the undeflected 
shaft position is noted. 

To measure the force at a particular current level, the magnet/rotor assembly is first 
degaussed using alternating current in the coils of the magnet, with peak amplitude of at 
least twice the highest current used. To avoid destructive vibrations while degaussing, the 
shaft is rigidly fixed relative to the magnets by using a temporary clamp. The magnets are 
then activated and the proximity probe outputs are read to determine the final shaft position 
relative to the magnets and also the absolute shaft deflection. 

3.3.3 Measurements Using One Magnet 

The apparatus described above was used to measure the force between a single magnet 
and the shaft for a variety of positions of the shaft with respect to the center of curvature of 
the magnet pole faces. Particular attention was paid to the forces when the shaft was given 
an eccentricity with respect to the axis of symmetry of the magnet. Such relative positions, 
which will be seen as undesirable, may nevertheless result from three causes: 
misalignment of the magnets during assembly (note that this is a strong argument for 
manufacture of magnets having poles attached to a continuous backing ring), from dynamic 
motion of the shaft, or from errors in biasing. 

Measurements were made at several levels of current in the magnet coils, and over a 
range of shaft positions within the clearance space. After assembly it was found that the 
magnets lacked a common center because of assembly tolerances. All the measurements 
therefore were conducted using one of the side magnets. At present the numerical 
calculation method described above has not been applied to the geometry of the new 
apparatus, so the results presented below are measurements only. A limited discussion of 
the trends of the normal forces in relation to the previous experiment as well as to the 
calculations that have been performed will be attempted, however. 

The results of these measurements are presented in a slightly different way from the 
results of the Section 2. 1 above, reflecting a different sequence of shaft repositioning from 
that used in the first series of experiments. Each group of symbols corresponds to a 
constant x position and therefore represents a traverse of the vertical direction (the normal 
direction in this case). By executing traverses of the normal direction it was possible to 
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plot the ratio of normal force to principal force as a function of the distance away from the 
axis of symmetry. 

In the plots below, the value of X, or x/c, listed in the legend gives the position of the 
rotor along the axis of symmetry of the magnet. The largest possible value is 1.0, but it 
was not possible to approach this value closely with the present design of the apparatus. A 
discussion of possible redesign to alleviate this difficulty, as well as to achieve some other 
goals, is presented in a later section. For the present, however, most of the values of X in 
this and subsequent plots are negative, and some extrapolation is necessary to visualize 
trends as the rotor approaches the magnet. Some of the values of y/c that are presented are 
smaller than - 1 .0, corresponding to a location so far from symmetry axis that it is outside 
the clearance space. While this would not be possible in an actual design, the steady state 
apparatus can accommodate such large displacement. 

Figure 3.3.2 shows the measured force in the principal direction and Figure 3.3.3 
shows the force in the normal direction at a current level of 1.0A (400 A t). The principal 
force is seen to be significantly larger at larger values of X, or rotor locations closer to the 
magnet. The force increases slightly as the rotor is moved away from the axis of 
symmetry, toward one of the poles. This trend was predicted by the original linear theory 
based on gap regions only, but has not been predicted by the nonlinear theory. The 
magnitude of the normal force increases strongly with an increase in distance away from 
the symmetry axis. Figure 3.3.4 shows the ratio of the normal force to the principal force. 
Linear fits to the data of each traverse were calculated. To avoid confusion only one fit is 
shown, but the slope of this line is equal to the average of the slopes of all the fits. 

Figures 3.3.5 through 3.3.16 show the corresponding results for other values of coil 
currents. Similar trends are observed in all the cases. In each of the plots of force ratio a 
linear fit is provided, and in each case the slope of the fit chosen is the same as the average 
of the slopes of all the individual fits. 

Over the ranges of position examined, the data indicate an approximately linear 
relationship between the normal eccentricity of the shaft and the ratio of normal to principal 
force. The constant of proportionality seems to be larger at lower currents, but for all cases 
examined its value is between 0. 14 and 0. 17. The nonlinear theory has predicted the 
existence of normal forces, but has not predicted such a large constant of proportionality 
for the ratio. 
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Figure 3.3. 1 1 Principal force from single magnet at 1.5 Ampere (600 A-turns). 





















92 


3.3.4 Measurements Using Two Magnets 

The result that the normal force seems to be approximately proportional to the 
product of principal force and normal eccentricity indicated that the normal forces would be 
even more significant in the case of a strongly biased, opposed magnet pair, where the 
principal forces are approximately balanced. To examine this possibility a new set of 
experiments was conducted using two opposed magnets at the same current. 

Figure 3.3.17 is a plot of the principal force exerted by the magnet pair as a 
function of the normal coordinate, at several values of the principal coordinate, along the 
symmetry axis. At small eccentricity, the force in the principal direction is near zero, as 
expected, since the shaft is equidistant from the two magnets. As the shaft is moved 
toward either of the magnets, there is a resultant force in the direction of that magnet. This 
force is a relatively weak function of the normal coordinate, but is seen to be largest in 
magnitude when both normal and principal eccentricities are large. 

The normal forces exerted under these conditions are shown in Figure 3.3. 18. The 
intercept of these force plots with the axes would nominally be at (0,0), and the precise 
cause of their displacement from that intercept is not yet clear, although its most likely 
cause seems to be uncertainty in the angular positions of the magnets around the clearance 
circle. This possibility is discussed further below, but it is felt that the magnitudes of the 
forces observed in this plot remain significant after considering this uncertainty, because all 
the normal forces do in fact change sign as anticipated, at some point. The magnitude of 
the normal force is the same order as the resultant principal force at small eccentricities and 
is significantly larger than the principal force at large eccentricities, as illustrated by the plot 
of Figure 3.3.19, which shows the ratio of normal to principal force. This occurs because 
the normal forces from the two magnets are additive, while the principal forces are of 
opposite sign. 

The displacement from zero of the intercepts of the normal forces with the axes in 
Figure 3.3. 19 and the locations of minima of the principal forces in Figure 3.3.17 have 
been considered in terms of possible angular uncertainty in magnet placement. The 
apparatus, further described in the earlier progress report, was assembled by positioning all 
of the magnet pole faces against a plastic mandrel and then clamping the magnet 
laminations in place. The mandrel was then carefully removed. Because of the method of 
assembly, it is felt that the uncertainty in radial position of each pole face is small, on the 
order of 1 % of the clearance. The largest uncertainty is that of the angular position of each 
magnet. 
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Ideally, each magnet should be 90° from each of its neighbors, but in fact, assembly 
tolerances may have resulted in errors of up to 3° from the nominal positions. This would 
result in an uncertainty in orientation of the force vector associated with each magnet. Such 
a 3° uncertainty would result in an increase or decrease in the normal force component of 
approximately 2 % of the principal force from each magnet of the pair. Using a principal 
force of 1 .0 at zero of the principal coordinate, a 3° uncertainty in principal force 
orientation would result in a change in the normal force of 0.06 nondimensional force 
units. This is an uncertainty of the same order as the displacement of intercepts of the 
measured normal force data from the (0,0) position in Figure 3.3.18. 

If this uncertainty in angular position is the cause of the displacement of the data, it 
would be expected that forces measured at other current levels would behave similarly. 
Additional measurements were made at two other current levels, and indeed the data all are 
displaced by approximately the same amount. These data are shown in Figures 3.3.20 
through 3.3.22, for MMF of 300 A t, and in Figures 3.3.23 through 3.3.25, for MMF of 
500 A t (each magnet). It is therefore believed that a small error in angular positioning is 
present. 

In interpreting the measurements of Figures 3.3.17 through 3.3.25, then, it must 
be remembered that this uncertainty of angular orientation may be playing a role. 
Nevertheless, the magnitudes of the normal forces are still significant. Those on the 
negative side of the plots are apparently increased by the error, but those on the positive 
side are apparently decreased. Therefore, it may be conservatively stated that the 
magnitudes of normal forces in a magnet pair that is perfectly aligned will be at least as 
large as those on the positive side of the plots presented. Under this interpretation the 
normal forces are of significant magnitude. In addition, the very fact that the system is 
shown to be highly sensitive to such angular errors should receive some emphasis. This 
factor must be considered in designing the actuator and in determining the level of 
robustness or the type of algorithms required for controlling the bearing. It is a strong 
argument, in fact, for the use of magnets made with a continuous outer ring, which would 
practically preclude this type of uncertainty. 

The coordinate coupling illustrated in the force ratio plots has serious implications 
for the design of magnetic bearings and magnetically supported flexible rotor systems 
because it is primarily dependent not on the control currents, but on the bias currents. 

Since it is widely believed that magnetic bearings must be biased rather strongly in order to 
provide a greater degree of linearity and to improve their stability characteristics, the 
coordinate coupling in these systems could be strong. In flexible rotor dynamics. 
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coordinate coupling in bearings is regarded as an undesirable characteristic because of the 
potential for excitation at multiples of the running speed as well as for self-excited whirl. 
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The type of coupling illustrated by these measurements would not tend to cause whirl, 
because the coupling coefficients have the same sign, unlike the case of a fluid film 
bearing, where the normal stiffness coefficients often have opposite signs. They would, 
however, tend to cause an excitation at multiples of the running speed. This possibility 
must be considered when designing magnetic bearings for flexible rotor applications, such 
as gas turbines and other turbomachinery. 
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4. CLOSURE 

4.1 Summary of Technical Findings 

Both numerical and experimental studies were done to determine the characteristics of 
the forces exerted on a stationary shaft by a magnetic bearing actuator. The numerical 
studies were based on finite element computations and were of three types: 

(1) Calculations based only on energy stored in air gaps. 

(2) Calculations including metal parts, with locally constant magnetic 
permeability. 

(3) Calculations including nonlinear magnetization function, with the possibility 
of saturation . 

Principal conclusions from the analytical/numerical studies are: 

(1) The distribution of saturation in the magnet core and the rotor influences both 
the principal attractive force and the normal force 

(2) The computed normal force is considerably smaller than that measured 
experimentally 

(3) The trend of nonlinear theory to predict larger normal forces in relation to 
principal forces is appropriate. 

Measurements of the force versus position of the shaft were made using two 
separate measurement rigs, one based on strain guage measurement of forces, the other 
based on deflections of a calibrated beam. All measurements were static, using steady 
currents and a nonrotating shaft. Principal conclusions from the experimental studies, 
taken in conjunction with the numerical studies, are: 

(4) The trends of the measured principal (y) forces agree with the predictions of 
the theory while the magnitudes of forces are somewhat smaller than those 
predicted. The y forces do not appear to decrease significantly as predicted 
by theory when the rotor is given an eccentricity in the x direction. 

(5) The measured forces in the x direction are significantly larger than those 
predicted by theory when the rotor has an x eccentricity. 

(6) Over the ranges of position examined, there is an approximately linear 
relationship between the normal eccentricity of the shaft and the ratio of 
normal to principal force. The constant of proportionality was not the same 
for all cases, but its value was consistently between 0. 14 and 0. 17. The 
nonlinear computations predicted the existence of normal forces, but did not 
predict such a large ratio. 
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(7) The type of coupling illustrated by these measurements probably would not 
tend to cause whirl, but they might tend to cause other kinds of 
nonsynchronous excitation This possibility must be considered when 
designing magnetic bearings for flexible rotor applications, such as gas 
turbines and other turbomachinery. 

Further work was conducted beyond the nominal work period of this proposal, in 
which simulations of 2DOF systems were performed subject to these force models. The 
results, attached as Appendix C, show that significant nonlinear behavior can occur, 
including multiple coexisting solutions, bifurcations in response as the stabilities of the 
respective solutions change, and self-similarity in stability boundaries. 

4.2 Documentation 

In the course of this project, one master's thesis and one Ph.D. dissertation were 
completed. Edward McCaul received the degree Master of Science after defending the 
thesis entitled "Measurement of Forces in a Magnetic Journal Bearing" [72], Harold Xia 
received the degree Doctor of Philosophy after defending the dissertation "Numerical 
Investigation of Suspension Force in a Magnetic Journal Bearing Actuator" [60], 

Five interim progress reports were filed with NASA as this work proceeded 

Presentations of results of the research performed under this grant were made at three 
technical meetings: 

1 . NASA Workshop on Aerospace Applications of Magnetic Suspension at Langley 
Research Center, September 1990 

2. ROMAG'91 Conference on Magnetic Bearings and Dry Gas Seals, in Alexandria, 
VA, March 1991 

3. ASME/STLE Joint Tribology Conference in St. Louis, Missouri, October 1991 
(also published in ASME Journal of Tribology) [68] 

In addition, related work on nonlinear dynamic simulation of magnetic bearing 
systems that makes direct use of the results obtained in this project have been presented at 
three technical conferences: 

1 . Third International Symposium on Magnetic Bearings, Alexandria, VA 1992 [69] 

2. NASA Second International Symposium on Magnetic Suspension Technology, 
Seattle, WA, August 1993 [73] 

3. ASME International Gas Turbine and Aeroengine Congress, The Hague, 1994 
(also published in ASME Journal of Engineering for Gas Turbines and Power) [74] 
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APPENDIX A - PROGRAM FOR FORCE CALCULATION USING 
AIR GAPS ONLY 

c$$$$$$$$$$$$$$$$$$$$$$$$$$$ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

$ 

c 

c GAPFOR SEPT. 1990 

C 

C THIS ALGORITHM CALCULATES THE MAGNETIC SUSPENSION FORCE FOR A 
C MAGNETIC JOURNAL BEARING ACTUATOR. THE CALCULATION REGION 
INVOLVES 

C THE AIR GAP BETWEEN THE MAGNET POLE-FACES AND THE ROTOR WHERE THE 

C MAGNETIC FIELD IS DOMINATED BY LAPLACE EQUATION. MAGNETIC FORCE 

C IS DETERMINED BY USING VIRTUAL WORK PRINCIPLE. 

C 

c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

$ 


c 

c INPUT 

c 

C LINE1-TITLE 

C LINE2- E: ECCENTRICITY OF SHAFT CENTER 

C C: CLEARENCE(IN. ) 

C R : JOURNAL RADIUS ( IN . ) 

C ANG1 , ANG2 : LOCATION OF THE MAGNET POLE 2 (DEGREES) 

C EBC: CODE USED FOR EXPANDING PART AS A PERCENTAGE OF THE 

C GAP REGION 

C (A) 0: 10% OF THE GAP REGION 

C (B) 1: 20% OF THE GAP REGION 

C CURRENT: CURRENT IN COILS (A) 

C NCOIL: # OF TURNS (EACH COIL) 

C LINE3- N: CIRCUMFERENTIAL ELEMENT DIVISIONS 

C M: RADIAL ELEMENT DIVISIONS 

C ALPHA: ANGLE BETWEEN POLES 

C NBON : SPECIFIED BOUNDARIES 

C PHISTAR1: ARBITRARILY ASSUMED B.C ALONG POLE 1 

C PHIR: BOUNDARY CONDITION ALONG ROTOR CURVATURE 

C MG_DETH : MAGNET DEPTH 

C 
C 

c OUTPUT 

c 

C PHI MAGNETIC POTENTIAL FIELD 

C BETA MAGNETIC FLUX DENSITY 

C GAMMA MAGNETIC FLUX 

C SIGMA MAGNETIC ENERGY 

C FORCE MAGNETIC SUSPENSION FORCE 

C 

c SUBROUTINES 

C 

C COORD GRID GENERATION 

C ELEMOD ELEMENT NODAL SPECIFICATION 

C ASSEML GLOBAL BANDED COEFFICENT MATRIX 

C MODMAT MODIFIED BANDED MATRIX 

C GAUSS2 GAUSS ELEMINATION METHOD FOR BANDED MATRIX 

C FLUXEG MAGNETIC FLUX U ENERGY 

C 
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$ 

c 

IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 

REAL *8 KAPPA1 , KAPPA2 , MMF , MG_DETH , MU 

REAL *8 BETA ( 1000 ) , AL (1000) , PHI (800) , PHIS (800) , XJJ (50),YJJ(50,50) 
REAL* 8 X(800) , Y ( 800 ) ,FJJ(50, 50) , RECMAT ( 800 , 100 ) , REMMAT (800, 100) 
REAL* 8 RX(300, 10) ,RY(300, 10) , RPHI ( 300 , 10 ) , XPLOT ( 500 , 4 ) 

REAL *8 YPLOT(500, 4) , PHIPLOT(500, 4) , BETAPLOT ( 1000 ) 

REAL* 8 ANGPLOT ( 1000 ) ,XBETA(1000) ,YBETA(1000) 

INTEGER EB, EBC, EM, ELEMT, B1 , B2 

INTEGER X_NUMBER , X_POINT , Y_POINT , Y_P 

INTEGER 11(1000) , JJ(1000) ,KK(1000> , IBOUND(2, 800) 

IN=1 

10=9 

C 

WRITE (IO, 5) 

5 FORMAT (1H1,10X, 1 * * * MAGNETIC FORCE CALCULATION — LINEAR MODEL 
* * * > ) 

READ (IN, 7) TITLE 
7 FORMAT ( 3 0A4) 

READ ( IN, 10 ) E , C , R , ANG1 , ANG2 , EBC , CURRENT , NCOIL 
10 FORMAT (5F10 . 4 , 14 , F10 . 4 , 16 ) 

WRITE (*,20)E,C,R, ANG1 , ANG2 

20 FORMAT (IX, ' E= ' , F5 . 3 , IX, ' C= ' , F5 . 3 , IX, ' R= 1 , F5 . 2 , IX, 

& ' ANG1= ' , F5 . 2 , IX , ' ANG2= ' ,F5.2) 

WRITE (*,25) EBC , CURRENT , NCOIL 

25 FORMAT ( IX, ' EBC= ' , 14 , IX, ' CURRENT= ' , F8 . 4 , IX, 'NUMBER OF COIL=',I6) 
READ ( IN, 30)N,M, ALPHA, NBON, PHISTAR1, PHIR, MG_DETH 
30 FORMAT (2110 , F10 . 4 , 110, 3F10 . 4) 

WRITE ( * , 35 ) M, N, ALPHA , NBON 

35 FORMAT (IX, ' M= ' , 14 , IX , ' N= ' ,14, IX, ' ALPHA= 1 ,F7.2,1X, ' NBON= ' ,14) 
WRITE ( * , 37 ) PHISTAR1 , PHIR , MG_DETH 
37 FORMAT (IX, ' PHISTAR1= ' ,F7.2,1X, ' PHIR= ' , F7 . 2 , 

& IX, 'MAGNET DEPTH=* ,F7 .4) 

CONVER=0 . 0254D0 
E=E*CONVER 
C=C*CONVER 
R=R*CONVER 

MG_DETH=MG_DETH * CONVER 
PI=4 . 0D0 *ATAN ( 1 . 0D0 ) 

MU=4 . *PI*1 . OE-7 

D_A= (ANG2-ANG1 ) * PI/ 180 .DO 

POL E_L = D_A * (R+C) 

WRITE (*, 55) POLE_L 
55 FORMAT ( IX , 1 POLE_L= ' , F15 . 9 ) 

IF ( EB . EQ . 0 ) GOTO 40 
EB=5 
NR=N/EB 
GOTO 50 
40 EB=10 

NR=N/EB 
50 Bl =NR* 2 + 1 
B2=N+B1 
N1=B2+NR*2 
Ml=M+l 

C TOTAL NODAL POINTS ( INCLUDING EXPANDED BOUNDARY AREA) 

NODE=Nl * (Ml) 
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N2=N1-1 

C TOTAL ELEMENTS 
ELEMT=2 *N2 *M 
WRITE ( * , 8 0 ) NODE , ELEMT 

^80 FORMAT ( IX, 'NODE=',H0, IX, 


1 ELEMT= 1 , HQ) 


c $$$$$$$$$ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 


DELTA=C / 1 0 .DO 
CHANGE=CONVER*l . OD-7 
DO 3010 1=1,8 
XJJ (I) =. 0D0 


3010 


C 


3020 

C3020 

3022 

C 


DO 3010 J=1 , 15 
YJJ ( I , J) = . 0D0 


DO 3020 1=1,1 


XJJ ( I ) = ( I-l ) * DELTA 
XJJ ( I) = . 0135D0*CONVER 




YJJ ( I, J) = . 021D0 *CONVER- ( J. 
WRITE ( * , 3022 ) XJJ (I) , YJJ (I 
FORMAT (2X, 1 XJJ= • , F12 . 8 , 2X, 


) * DELTA 
) 

YJJ= ' , F12 . 8 ) 


C LOOP ON FORCE CALCULATION 


DO 1050 1=1,1 
DO 1050 J=l,l 





$$$$$$$$$ 


XJ=XJJ ( I ) -CHANGE 
YJ =YJJ (I, J) 
A1=ANG1+ ALPHA 
A2 = ANG2 + ALPHA 


CALL COORD ( Bl, B2, HI, EB,N,M,C,R,A1,A2 - PI.XJ, YJ, X, Y, 


RANG1 , RANG2 , 


CALL 

C 

C 

C CHECK ORDER OF VERTICES 
WRITE(IO, 8000) 

noT!2; '^ E L K °? DER 0F VERTICES WITHIN 


ELEMOD ( M , N2 , I I , J j , KK , ELEMT ) 


WITHIN ELEMENTS 


8000 


82 


DO 82 ITEST=1, ELEMT — — wxxhxjm ELEMENTS') 

WRITE ( IO, 84 ) II ( ITEST) , JJ(ITEST) KK I TTFc-t' \ 

' ' 11 ( ITEST) > ' Y ( JJ ( ITEST) ) , Y(KK (ITEST) ) 
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84 FORMAT (IX, (316, IX) ) 
c 8 5 bandwidth' of' global' RECTANGULAR MATRIX 
NWID=Ml+2 
WRITE (* , 135)NWID 

13 5 FORMAT ( IX, 'NWID=' ,16) 

c 

c 


call 

Sc 


c— 


c 

c DO 110 K=1 , NODE NWID).PHI(K) 

CLIO WRITE!*. 120) CRECMATtK.KJI 

120 FORMAT ( IX ,7F8.2,2X,F8. 

Tall 

& c IBKK,IBJJ,M,N,EB) ^ 


C " " _ CALL GAUSS2 ( REMMAT , PHI , NODE ,NWID)^^ 

C WRITE (IO, 140) Ml, N1 

140 FORMAT(3X, 15, ' , ' .I 5 ) 

MNODE=NODE-M 

r* DO 150 IPOT = l , Ml 

r DO 160 JPOT=IPOT , MNODE , Ml 

c X ( JPOT) =X ( JPOT) /CONVER 

C C 1S„ .RHI.OPOT, 

C150 MNODE=MNODE+ 1 3) 

170 FORMAT (3X, 2 (F7 . 3 , , ) . • J 

c 


TT TT KK RANG1,RANG2,N,AL, PHI, X,Y, GAMMA, 




gammastari=gamma 
C WRITE ( * , 220 ) GAMMASTARl ^ 

220 FORMAT (IX, 'GAMMASTARl- ,F12. 

PHI STAR2 = - PHI STARl 

A1=ANG1 

A2 =ANG2 

C 

C-- 


, „ t vt y v RANG 1 , RANG 2 , 

CALL COORD(Bl,B2,Nl,EB,N,M,C,R,Al,A2,P 

& RX,RY) 
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c 

c 

C 

CALL ELEMOD (M,N2, II, JJ, KK , ELEMT ) 

c 

C 

c 

CALL ASSEML (M, N, NWID, ELEMT, X, Y, EB , NODE , II , JJ , KK, RECMAT , PHI , AL , 

& PHISTAR2, PHIR, PHIS, IBKK, IBJJ, IBOUND) 

C 

C 

c 

CALL MODMAT ( NWI D , NBON , PHI , NODE , RECMAT, REMMAT, PHIS, IBOUND, 

& IBKK, IBJJ, M,N,EB) 

C 

C 

C 

CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 

C 

c 

c 

CALL FLUXEG ( EB , ELEMT , II , JJ , KK , RANG1 , RANG 2 , N , AL , PHI , X , Y , GAMMA , 

$ SIGMA , BETA , M, R, C , NODE , MG_DETH , ANGPLOT , XBETA , YBETA) 
c 

C 

C DO 230 I=EM, 1 , -M2 

C230 WRITE ( * , 240 ) (BETA(J) ,J=I,I+M2-1) 

C 

c 

C CALCULATE PARAMETERS TO DECIDE BOUNDARY CONDITION 
C 

c 

GAMMA STAR 2 =GAMMA 
C WRITE (*, 250 )GAMMASTAR2 

250 FORMAT ( IX, ' GAMMASTAR2= ' , F12 . 4 ) 

KAPPA1=GAMMASTAR1 / PHISTAR1 
KAPPA2=GAMMASTAR2 / PHISTAR2 
MMF=CURRENT* FLOAT (NCOIL) 

PHIP1=MMF/ (1 . -KAPPA1/KAPPA2 ) 

PHIP2=PHIP1-MMF 

C 

c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

$ 

C CALCULATION OF MAGNETIC FLUX AND ENERGY FOR BACKWARD PERTURBATION 
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 
$ 

C 

XJ=XJJ ( I ) -CHANGE 
YJ=YJJ ( I , J) 

A1=ANG1+ ALPHA 
A2=ANG2+ ALPHA 
C 

c 
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NK=1 

DO 650 IPLOT=l,N2 

DO 660 JPLOT=l,M 

XPLOT (NK, 1 ) =RX ( IPLOT+1 , JPLOT) 

XPLOT (NK, 2 ) =RX ( I PLOT, JPLOT) 

XPLOT ( NK , 3 ) =RX ( I PLOT , JPLOT+ 1 ) 

XPLOT (NK, 4 ) =RX ( IPLOT+1 , JPLOT+1 ) 

YPLOT (NK , 1 ) =RY (IPLOT+1 , JPLOT) 

YPLOT ( NK , 2 ) =RY ( I PLOT , JPLOT ) 

YPLOT (NK, 3 ) =RY ( IPLOT , JPLOT+1 ) 

YPLOT (NK, 4) =RY( IPLOT+1, JPLOT+1) 

PHIPLOT (NK, 1 ) =RPHI ( IPLOT+1 , JPLOT) 

PHI PLOT (NK, 2 ) =RPHI ( IPLOT , JPLOT) 

PHIPLOT (NK, 3 ) =RPHI ( IPLOT , JPLOT+1 ) 

PHIPLOT (NK, 4 ) =RPHI ( IPLOT+1 , JPLOT+1 ) 

660 NK=NK+1 
650 CONTINUE 

NUMELE=ELEMT / 2 
C WRITE (IO, 670) NUMELE 

670 FORMAT ( 3X, 15 ) 

C DO 680 I PLOT=l, NUMELE 

C XPLOT (IPLOT!icON) =XPLOT (IPLOT , JCON) /CONVER 

C 6 8 5 YPLOT ( I PLOT , JCON ) =YPLOT (IPLOT, JCON ) /CONVER 

C DATA OUTPUT FOR CONTOUR GRAPHICS 

C$ $$$$$$$$$$$$$ $$$$$$$$$$$$$$$$$$$$$$ 

C 

c WRITE (IO, 8010) 

C8010 FORMAT ( 2X , 1 CONTOUR PLOT') 

c WRITE (IO, 690) (XPLOT ( IPLOT , JPLOT) , JPLOT-1,4 , 

Z. & (YPLOT (IPLOT, JPLOT) , JPLOT=l , 4) 

C680 WRITE (10,700) ( PHIPLOT ( IPLOT, JPLOT) , JPLOT=l , 4 ) 

690 FORMAT ( 3X, 7 (F7 . 3 , ' , ' ) , F7 . 3 ) 

700 FORMAT ( 3X, 3 (F10 . 4 , ' , ' ) , F10 . 4) 

Al=ANGl 

A2=ANG2 

CALL, C00RD(B1,B2,N1,EB,M,M,C,R,A1,A2,PI,XJ,YJ,X,Y,RANG1,RANG2, 

& RX f RY) 

C 

CALL ELEMOD (M, N2 , II , JJ, KK, ELEMT) 

C 

CALL ASSEML(M,N, NWID, ELEMT, X,Y,EB, NODE II, J J , KK , RECMAT , PHI, AL, 
& PHIP2 , PHIR, PHIS, IBKK, IBJJ, IBOUND)^ 
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CALL MODMAT { NWID , NBON , PHI , NODE , RECMAT , REMMAT , PHIS , IBOUND, 
& IBKK , IBJJ ; M,N, EB) 



C 

c 

CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 

c 

C 




CALL FLUXEG ( EB , ELEMT f I I , JJ , KK , RANGl , RANG 2 , N , AL , PHI , X , Y , GAMMA, 
$ SIGMA , BETA , M , R , C , NODE , MG_DETH , ANGPLOT , XBETA , YBETA ) 



C 

GAMMA_2 A=G AMMA 
S I GMA__ 2 A= SIGMA- 

WRITE ( * , 310) G AMMA_2 A , SIGMA_2A 
310 FORMAT (IX # 1 GAMMA_2A= ' , FI 2 . 4 , IX, ' SIGMA_2A= ' , F12 . 4 ) 

C 

C CALCULATION OF BOUNDARY CONDITIONS FOR FORWARD PERTURBATION 

C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

C 

XJ-XJJ ( I ) +CHANGE 
YJ=YJJ ( I , J ) 

A1=ANG1+ ALPHA 
A2 =ANG2 +ALPHA 
C 




CALL COORD ( B1 , B2 , N1 # EB , N , M , C , R , Al , A2 , PI , XJ , YJ , X , Y , RANGl , RANG2 , 
& RX , RY ) 


C 

C 

CALL ELEMOD (M, N2 , II ,JJ,KK, ELEMT) 

C 

c 

c 


CALL ASSEML(M,N,NWID, ELEMT, X, Y, EB , NODE, II , JJ , KK , RECMAT , PHI, AL, 
Sc PHISTARl , PHIR, PHIS, IBKK, IBJJ, IBOUND) 


C 



CALL MODMAT { NWID , NBON , PHI , NODE , RECMAT , REMMAT , PH I S , I BOUND , 
Sc IBKK, IBJJ , M , N , EB) 



C 

c 

CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 

C 

C 



nnio noio tooo 0010 100 
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CALL FLUXEG ( EB , ELEMT , I I , J J , KK , RANG1 , RANG 2 , N , AL , PHI , X , Y , GAMMA , 
$ SIGMA , B ETA , M, R, C, NODE , MG_DETH , ANGPLOT , XBETA , YBETA ) 
c 

c 

M2=M*2 

EM=ELEMT-M2+1 
GAMMAS TAR 1 =GAMMA 
C WRITE (*, 410) GAMMASTAR1 

410 FORMAT (IX, 1 GAMMASTAR1 = ' ,F12.4) 

PHISTAR2=-PHISTAR1 

Al=ANGl 

A2=ANG2 


CALL COORD (B1 , B2 , N1 , EB , N, M, C , R, A1 , A2 , PI , XJ , YJ , X, Y , RANG1 , RANG2 , 
& RX, RY) 


CALL ELEMOD(M,N2, II, JJ,KK, ELEMT) 


CALL ASSEML (M, N, NWID, ELEMT, X, Y, EB, NODE, II , JJ, KK, RECMAT , PHI , AL , 
& PHISTAR2, PHIR, PHIS, IBKK, IBJJ, IBOUND) 


CALL MODMAT (NWID, NBON, PHI , NODE , RECMAT, REMMAT, PHIS, IBOUND, 
& IBKK, IBJJ, M,N,EB) 


CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 

C 

C 

C 


CALL FLUXEG ( EB , ELEMT , I I , JJ , KK , RANG1 , RANG2 , N, AL , PHI , X , Y , GAMMA , 
$ SIGMA, BETA, M, R, C , NODE , MG_DETH , ANGPLOT, XBETA, YBETA) 
c 

C 

GAMMASTAR2 =GAMMA 
C WRITE (*, 420) GAMMASTAR2 

420 FORMAT (IX, ' GAMMASTAR2 = ' ,F12.4) 

KAPPA1=GAMMASTAR1 / PHISTAR1 
KAP PA2 =GAMMASTAR2 / PHI STAR2 
MMF=CURRENT* FLOAT (NCOIL) 

PHIP1=MMF / ( 1 . -KAPPA1/KAPPA2 ) 

PHIP2=PHIP1-MMF 

C WRITE (*,430) KAPPA1 , KAPPA2 , FP1 , FP2 



no oin inon non noin inoo ooin ion 
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430 FORMAT (IX, 'KAPPA1=' , F12.4,1X, 'KAPPA2=' ,F12.4,1X, 

& ' FP1= ' , F12 . 4 , IX, ' FP2= 1 , F12 . 4 ) 

C 

c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

C CALCULATION OF MAGNETIC FLUX AND ENERGY FOR FORWARD PERTURBATION 
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 

c 

XJ=XJJ ( I ) +CHANGE 
YJ=YJJ ( I , J) 

Al=ANGl+ ALPHA 
A2=ANG2+ALPHA 


CALL COORD (B1 , B2 ,N1 , EB, N,M, C, R, A1 , A2, PI, XJ, YJ, X, Y, RANG1 , RANG2 , 
& RX , RY ) 


CALL ELEMOD (M,N2, II, JJ, KK , ELEMT ) 


CALL ASSEML (M, N, NWID, ELEMT , X , Y, EB , NODE, II , JJ, KK, RECMAT , PHI, AL, 
& PHIP1 , PHIR, PHIS, IBKK, IBJJ, IBOUND) 


CALL MODMAT (NWID, NBON, PHI , NODE, RECMAT, REMMAT, PHIS, IBOUND, 
& IBKK, IBJJ, M,N, EB) 


CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 


CALL FLUXEG ( EB , ELEMT , I I , JJ , KK , RANG1 , RANG 2 , N , AL , PHI , X , Y , GAMMA , 
$ SIGMA , BETA , M, R, C , NODE , MG_DETH , ANGPLOT , XBETA, YBETA ) 


GAMMA_1B= GAMMA 
SIGMA_1B=SIGMA 

WRITE (*, 330) GAMMA_1B, SIGMA_1B 
0 FORMAT (IX, ' GAMMA_1B= 1 ,F12.4,1X, ' SIGMA_1B= 1 , F12 . 3 ) 
A1=ANG1 
A2=ANG2 


CALL COORD ( B1 , B2 , N1 , EB , N, M , C , R , A1 , A2 , PI , XJ, YJ, X, Y, RANG1 , RANG2 , 
& RX , RY ) 



ooio inon on 


All 


c 


CALL ELEMOD (M,N2 , II , JJ, KK, ELEMT) 


CALL ASSEML ( M , N , NWID , ELEMT , X , Y , EB , NODE , I I , J J , KK , RECMAT , PHI , AL , 
PHIP2 , PHIR, PHIS, IBKK, IBJJ, IBOUND) 


C — 
C 

c 

c — 
c 

C-- 


CALL MODMAT ( NWID , NBON , PHI , NODE , RECMAT , REMMAT, PHIS, IBOUND, 
& IBKK, IBJJ, M,N, EB) 


CALL GAUSS2 (REMMAT, PHI, NODE, NWID) 


CALL FLUXEG (EB, ELEMT, II , JJ , KK , RANG1 , RANG2 , N, AL, PHI , X, Y, GAMMA, 
$ S IGMA , BETA , M, R, C , NODE , MG_DETH , ANGPLOT , XBETA , YBETA ) 


C 

GAMMA_2 B =GAMMA 
SIGMA_2B= SIGMA 

WRITE ( * , 3 4 0 ) GAMMA_2 B , S IGMA_2 B 

340 FORMAT ( IX, 1 GAMMA_2B= ' ,F12.4,1X, 'SIGMA_2B=' ,F12.3) 

C 

C FORCE CALCULATION 
C 

DELTA_SIGMA= (SIGMA_1B+SIGMA_2B) - ( SIGMA_1A+SIGMA_2A) 

FORCE=DELTA_SIGMA/CHANGE 

FORCE=MU*FORCE/4 . 

F J J ( I , J ) =FORCE 

WRITE ( * , 1070 ) XJJ ( I ) , YJJ (I , J) ,FJJ(I,J) 

1050 CONTINUE 

DO 1060 1=1,1 
XJJ(I)=XJJ(I)/C 
DO 1060 J=1 , 15 
YJJ (I, J) =YJJ (I, J) /C 
WRITE(IO, 8050) 

8050 FORMAT (IX, 4X, ' X ' , 12X, ' Y 1 , 17X, 'FORCE' ) 

1060 WRITE (IO, 1070) XJJ (I) , YJJ (I, J) , FJJ (I, J) 

1070 FORMAT ( IX, F12 . 8 , IX , F12 . 8 , IX, F12 . 3 ) 

STOP 

END 


C 

£************ 


*★★* + ****** + 


********************************************** 


SUBROUTINE COORD (Bl , B2 , N1 , EB , N, M, C , R, ANG1 , ANG2 , PI , XJ , YJ , X , Y , 
&RANG1 , RANG2 , RX , RY ) 
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Q********************************************************************** 


IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 

REAL* 8 X(800), Y(800),RY (300,10), H(300),RX (300, 10) 

INTEGER EB,B1,B2 
10=9 

RANG1=ANG1*PI/180 . 

RANG2=ANG2*PI/180 . 

RANG= (ANG2-ANG1) /FLOAT (EB) *PI/180. 

THETA1=RANG1-RANG 

THETA2 =RANG2 +RANG 

THETA=THETA1 

DO 20 1=1, Nl 

IF ( I . NE . 1 ) GO TO 30 

H ( I ) =C-XJ*COS (THETA) -YJ*SIN (THETA) 

GO TO 35 

30 IF ( I . LE . Bl ) GO TO 40 

IF ( I . GT . B2 ) GO TO 40 
DELTA_TH= ( RANG2 -RANG1 ) /FLOAT (N) 

GO TO 50 

40 DELTA_TH= ( RANG2 -RANG1 ) /FLOAT (N) *0 . 5 

50 THETA=THETA+DELTA_TH 

H ( I ) =C-XJ*COS (THETA) -YJ*SIN (THETA) 

35 CC=COS (THETA) 

DELTA_H=H ( I ) /FLOAT (M) 

DO 80 J=1 , M+l 

RX ( I , J ) = (R+C-H ( I ) +DELTA_H* ( J-l) ) *COS (THETA) 

RY ( I , J ) = (R+C-H ( I ) +DELTA_H* (J-l) ) *SIN (THETA) 

80 CONTINUE 
20 CONTINUE 
K=1 

DO 120 1=1, Nl 

DO 130 J=1 , M+l 

X(J+ (K-l) * (M+l) ) =RX ( I , J) 

130 Y(J+ (K-l) * (M+l) ) =RY(I, J) 

120 K=K+1 

RETURN 
END 

***★★★★**★★★********★★* + *********★★★★★★****★**★*■*■***★****★*•*:*********★ 

SUBROUTINE ELEMOD (M, N2 , II, JJ,KK, ELEMT) 
£******************************+*************************************** 


INTEGER 11(1000) , JJ(1000) ,KK(1000) , ELEMT, ENDE 

INTEGER STARTE, T, END 

10=9 

M1=M+1 

STARTE=1 

ENDE=STARTE+ ( M- 1 ) * 2 
K=1 

END=ELEMT-1 
20 CONTINUE 
T=0 

DO 30 I=STARTE , ENDE, 2 
II (I) =K+T 



A13 


JJ(I) =n (I) +M1 + 1 
KK ( I ) =JJ ( I ) -1 
II (1+1) =11 (I) 
JJ(I+1)=II(I+1)+1 
KK ( I + 1 ) = JJ (I) 

30 T=T+1 

IF ( ENDE . EQ . END) GO TO 40 
STARTE=ENDE+2 
ENDE=ENDE+M*2 
K=K+M+1 
GO TO 20 
40 CONTINUE 
RETURN 
END 


C**************** ****************************************************** 

SUBROUTINE ASSEML (M, N, NWID, ELEMT , X, Y, EB, NODE, II , JJ, KK, RECMAT, PHI , 
f, AL , PHISTARl , PHIR, PHIS , IBKK, IBJJ , IBOUND) 

***************************************************************** 


C 


10 

20 


30 


40 


C 

41 

Y(I) 


IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 

REAL* 8 X(800) , Y (800) ,B(3) ,C(3) , AL (1000) , PHI (800) 
REAL* 8 RECMAT (800 ,100) , BDBMAT (1000 ,3,3) 

REAL* 8 PHISTARl, PHIR, PHIS (800) 

INTEGER IBOUND (2 , 800) , ELEMT, EB 

INTEGER 11(1000) , JJ ( 1 000 ) , KK ( 1000 ) , INELE (3,1000) 
10=9 

DO 10 1=1, NODE 
DO 10 J=1 , NWID 
RECMAT (I , J) =0 . 0 
DO 20 1=1, NODE 
PHI (I) =0.0 
M1=M+1 

M2= (N/EB*2+1) *M1 
M3 =NODE-N/EB* 2 *M1 


IBKK=1 

DO 30 I=M2 , M3 , Ml 


PHIS ( I ) =PHISTARl 
IBOUND (1, IBKK) =1 
IBKK=IBKK+1 
IBKK=IBKK-1 
IBJJ=1 

DO 40 1=1, NODE, Ml 
PHIS ( I ) =PHIR 
IBOUND (2, IBJJ) =1 
IBJJ=IBJJ+1 
IBJJ=IBJJ-1 
DO 50 NN=1, ELEMT 
I=II (NN) 

J=JJ (NN) 

K--KK (NN) 

WRITE (IO, 41) X (I) ,X(J) , X (K) ,Y(I) ,Y(J) ,Y(K) 


FORMAT ( IX, ( 6F8 . 6 , IX) ) 

AL(NN) =0 .5* (X ( I ) *Y(J) +X(J) *Y(K)+X(K) *Y(I) -X(I) Y(K)-X(K) Y(J) 


&*X(J) ) 
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C WRITE (10, 42)AL(NN) 

42 FORMAT ( IX, F12 . 8 ) 

B(1)=Y(J) -Y(K) 

B (2 ) =Y (K) -Y ( I ) 

B ( 3 ) =Y ( I ) -Y ( J) 

C ( 1 ) =X ( K ) -X(J) 

C ( 2 ) =X ( I ) -X ( K) 

C ( 3 ) =X ( J) -X(I> 

INELE ( 1 , NN) =11 (NN) 

INELE ( 2 , NN) =JJ(NN) 

INELE ( 3 , NN ) =KK ( NN ) 

DO 60 NR=1 , 3 
DO 60 NC=1 , 3 

60 BDBMAT (NN, NR, NC ) = ( B (NR) *B (NC ) +C (NR) *C (NC) ) / ( 4 . 0*AL (NN) ) 

50 CONTINUE 

DO 70 NN=1 , ELEMT 
DO 80 IN=1 , 3 
DO 80 JN=1 , 3 
ID=INELE ( IN, NN) 

JD=INELE ( JN, NN) 

IF ( JD. LT. ID) GO TO 85 

RECMAT ( ID, JD+l-ID) =RECMAT ( ID, JD+l-ID) +3DBMAT (NN, IN, JN) 

GOTO 80 

85 DO 100 NEWB=1,IBKK 

I F ( JD . NE . I BOUND ( 1 , NEWB ) ) GOTO 100 

PHI (ID) = PHI (ID) -BDBMAT (NN, IN, JN) *PHIS ( JD) 

100 CONTINUE 
80 CONTINUE 
70 CONTINUE 
RETURN 
END 

********************************************************************** 

SUBROUTINE MODMAT { NWID , NBON , PHI , NODE , RECMAT , REMMAT , PHIS , 

Sc I BOUND , IBKK, IBJJ,M,N,EB) 

£********************************************************************** 

* 

C 

IMPLICIT DOUBLE PRECISION (A-H, O-Z ) 

REAL* 8 PHI ( 800 ) , RECMAT ( 800 , 100 ) ,PHIS(800) , REMMAT ( 800 , 100 ) 

INTEGER MM ( 10 ) ,NEND(10) , IBOUND ( 2 , 800 ) , IIBOUND (2 , 800 ) 

INTEGER BON1 , BON2 , EB 

10=9 

Ml=M+l 

M2= (N/EB*2+1) *M1 
M3=NUMN-N/EB*2 *Ml 
DO 40 1=1, NODE 
DO 40 J=1 , NWID 
40 REMMAT ( I , J ) =RECMAT ( I , J ) 

NEND ( 1 ) = I BOUND ( 1 , IBKK) 

NEND ( 2 ) =IBOUND ( 2 , IBJJ ) 

MM ( 1 ) =M2 
MM (2 ) =1 

DO 100 NN=1 , NBON 
BONl=MM (NN) 

DO 110 1=1, NODE 

IF (I .GT.NEND(NN) )GO TO 110 
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IF ( I. EQ. BONl) GO TO 130 
SUM=0 . 0 
BON2 =BONl 
JJ=NWID+I-1 
DO 140 J=I,JJ 

IF (BON2 . GT . NEND (NN) ) GO TO 140 
I F ( J . NE . BON2 ) GO TO 140 
SUM=SUM+REMMAT ( I , J+l-I ) *PHIS(J) 

BON2 =BON2 +M1 
REMMAT ( I , J+l-I) =0 . 0 
140 CONTINUE 

PHI (I) =PHI (I) -SUM 
GO TO 110 
130 CONTINUE 

DO 170 J=1 , NWID 
IF ( J . EQ . 1 ) GO TO 170 
REMMAT ( I, J) =0.0 
170 CONTINUE 

PHI ( I ) =REMMAT ( I , 1) *PHIS (I) 

BONl=BONl+Ml 
110 CONTINUE 
100 CONTINUE 
RETURN 
END 

******************************************** 


SUBROUTINE GAUSS2 ( E, V, NNODES , NWID) 
£***************★*******★*★****************** 
* 

c 

IMPLICIT DOUBLE PRECISION (A-H, O-Z ) 

REAL *8 E(800, 100) , V ( 800 ) , U ( 800 ) 

C 

C TRIANGULARIZE 

C 

NWIDM=NWID-1 
ISTOP=NNODES-l 
ISTOP=NNODES-l 
DO 1000 ID=1 , ISTOP 
DO 1000 JD=1 , NWIDM 
QUO=E ( ID, JD+1 ) /E ( ID, 1 ) 

V ( ID+JD) =V ( ID+JD) -QUO*V ( ID) 

KSTOP=NWID- JD 
DO 1000 KD=1 , KSTOP 

E (ID+JD, KD) =E (ID+JD, KD) -QUO*E ( ID , KD+ JD) 
1000 CONTINUE 
C 

C BACK SUBSTITUTE 

C 

U ( NNODES ) =V{ NNODES) /E (NNODES , 1 ) 

DO 3000 ID=2, NNODES 
SUM=0 . 

IN=NNODES+l-ID 
DO 2500 JD=2 , NWID 
SUM=SUM+E ( IN, JD) *U(IN+JD-1) 

2500 CONTINUE 

U ( IN) = ( V ( IN) -SUM) /E ( IN, 1 ) 
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3000 CONTINUE 

DO 4000 ID=1 , NNODES 
V ( ID ) =U ( ID) 

4000 CONTINUE 
RETURN 
END 


************** 


r***************************** 


************************** 


SUBROUTINE FLUXEG { EB , ELEMT , II , JJ, KK, RANG1 , RANG2 , N, AL, PHI1 , X, Y, 
$GAMMA, SIGMA, BETA, M, RR, CL, NODE , MG_DETH, ANGPLOT , XBETA, YBETA) 

c , ********************************************************************** 


C 

IMPLICIT DOUBLE PRECISION ( A-H , O-Z ) 

REAL* 8 AL (1000 ) , PHI1 (800),X(800),Y(800) ,MG_DETH, ANGPLOT ( 1000 ) 
REAL* 8 BETA ( 1000 ) ,BETAX(1000) ,BETAY(1000) ,B(800) ,C(800) 

REAL*8 XBETA ( 1000 ), YBETA ( 1000 ) 

REAL * 8 RANG1 , RANG2 , GAMMA , SIGMA , RR , CL , DELTA_TH , DELTA_L 
INTEGER ELEMT, EB, II (1000) , JJ (1000) 

INTEGER KK(1000) 

10=9 

DO 20 1=1, ELEMT 
BETAX ( I ) =0 . 0 
20 BETAY ( I ) =0 . 0 
C 

C CALCULATE ELEMENT CENTER 
C 

DO 25 ICENT=1 , ELEMT, 2 
IE=II (ICENT) 

JE=JJ ( ICENT) 

KE=KK( ICENT) 

XMID= ( X ( IE ) +X ( KE ) )/ 2. 

YMID= ( Y ( IE ) +Y ( KE ) ) / 2. 

XBETA ( ICENT) = (XMID+0 . 5*X ( JE) ) /1.5 
YBETA ( ICENT) = ( YMID+0 . 5*Y ( JE) )/1.5 
IE=II ( ICENT+1 ) 

JE=JJ ( ICENT+1 ) 

KE=KK( ICENT+1) 

XMID= (X ( IE) +X (KE) ) 12 . 

YMID= (Y (IE) +Y (KE) )/ 2. 

XBETA ( ICENT+1 )= (XMID+0 . 5*X ( JE) ) /I . 5 
YBETA ( ICENT + 1 ) = (YMID+0 . 5 * Y ( JE ) ) / 1 . 5 
25 CONTINUE 
C 

C CALCULATE FLUX DENSITY 
C 

DO 30 NN=1, ELEMT 
I=II (NN) 

J=JJ (NN) 

K=KK (NN) 

B ( I ) =Y ( J ) -Y(K) 

B ( J ) =Y ( K) -Y ( I ) 

B(K)=Y(I)-Y(J) 

C(I) =X(K) -X(J) 

C ( J ) =X ( I ) -X (K) 

C ( K ) =X ( J ) -X ( I ) 

BETAX (NN) =- (B (I) *PHI1 (I) +B( J) *PHI1 (J) +B(K) *PHI1 (K) ) / (2 . *AL (NN) ) 



BETAY(NN) =-(C(I) *PHIl (I ) +C ( J) *PHIl (J)+C(K) *PHI1 (K) ) / (2 . *AL (NN) ) 
C WRITE (*, 100 ) BETAX (NN) ,BETAY(NN) 

100 FORMAT (2X, ’ BX= ' , F15 . 4 , IX, 1 BY= ' , F15 . 4 ) 

ANGPLOT ( NN ) =ATAN ( BETAY ( NN ) / BETAX ( NN ) ) 

30 BETA (NN) =SQRT (BETAX (NN) **2+BETAY (NN) **2) 

C WRITE (*,110) 

110 FORMAT (2X) 

DELTA_TH= ( RANG2 -RANG1 ) /FLOAT (N) *0.5 
DELTA_L=2 . * (RR+CL) *SIN (DELTA_TH) 

C DELTA_L=0 . 1 

C WRITE (*, 38) DELTA_TH,DELTA_L 

38 FORMAT ( IX, ' DELTA_TH= ' , F12 . 9 , IX, ' DELTA_L= 1 , F12 . 9 ) 

GAMMA=0 . 0 
M2=M*2 

NSTART=M2 * ( N / EB * 2 + 1 ) 

NEND=ELEMT-N/EB*2 *M2 
DO 40 I =NSTART , NEND , M2 
4 0 GAMMA=GAMMA+ BETA ( I ) 

GAMMA = DELTA_L * GAMMA * MG_DETH 
C GAMMA=GAMMA*MG_DETH 

SIGMA=0 . 0 
DO 50 1=1 , ELEMT 

50 SIGMA=SIGMA+ (BETAX ( I ) * *2+BETAY ( I ) * *2 ) *AL ( I ) *MG_DETH 
RETURN 
END 


C 

£****************** 


************************************ 


************* 


SUBROUTINE GAUSS 1 < E , V, NNODES , NWID) 

C***************************************************** 

C 


************** 


IMPLICIT DOUBLE PRECISION (A-H , O-Z ) 
REAL * 8 E(200,200) ,V(800) , U ( 800 ) 


C 

C TRIANGULARIZE 

C 

NWIDM=NWID-1 
DO 1000 ID=1, NNODES 
DO 1000 JD=1 , NWIDM 
QUO=E ( ID+JD, ID) /E ( ID, ID) 

IDP=ID+NWIDM 

V ( ID+JD) =V ( ID+JD) -QUO*V ( ID) 

DO 1000 KD=ID, IDP 

E ( ID+JD, KD) =E ( ID+JD, KD) -QUO*E ( ID, KD) 
1000 CONTINUE 
C 

C BACK SUBSTITUTE 

C 

U ( NNODES ) =V ( NNODES ) / E ( NNODES , NNODES ) 
DO 2000 ID=2 , NWID 
SUM=0 . 

IN=NNODES+l-ID 

NSTOP=NNODES 

INP=IN+1 

DO 1500 JD=INP , NSTOP 
SUM=SUM+E ( IN , JD) *U(JD) 

1500 CONTINUE 

U ( IN) = (V ( IN) -SUM) /E(IN, IN) 

2000 CONTINUE 



NWIDP=NWID+1 

DO 3000 ID=NWIDP , NNODES 

SUM=0 . 

IN=NNODES+l - ID 
NSTOP=IN+NWID- 1 
INP=IN+1 

DO 2500 JD=INP, NSTOP 
SUM=SUM+E ( IN, JD) *U ( JD) 

2500 CONTINUE 

U ( IN) = {V ( IN) -SUM) /E(IN, IN) 
3000 CONTINUE 

DO 4000 ID=1, NNODES 
V ( ID ) =U ( ID) 

4000 CONTINUE 
RETURN 
END 
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APPENDIX B - COMPUTATIONAL METHODS INCLUDING METAL 

B.l Modelling of Magnetization Curve 

If hysteresis and anisotropic effects are neglected, the magnetization function, or B 
vs. H function, for a ferromagnetic material can be represented as a continuous, single- 
valued nonlinear function. It is difficult, however, to find a single analytical expresston 
t at accurately represents the function over the entire useful range. A procedure using 
cubic sphne fits recommended by Silvester et al.fBI] is used in the present work to model 
the magnetization characteristic of silicon sheet steel. The graphical data for the 
characteristic were taken from Smith [B2], 

The procedure for modelling the experimentally determined relationship between B 
and H is as follows: 

First, convert the data from a permeability representation (B vs. H) to a reluctivity v 

vs. B2. the square of the flux density. The two methods of representing the data are shown 
in Figure B1 and Figure B2. 


Silicon sheet steel 
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Figure B2. Reluctivity vs. square of flux density 


Next, evenly divide B 2 into n subintervals. Then within an interval between end 
points i 1 and i2 the value of v and the slope K of the v-B 2 curve can be calculated by the 
cubic interpolation formulae 


v(x) = (2x 3 -3x 2 +l) Vn+(-2x 3 +3x 2 ) Vj2+(x 3 -2x 2 +x) Kn+(x 3 -x 2 ) Kj 2 (1.2) 

k(x) = -^- = — -J— -{(6x 2 -6x) Vji+(-6x 2 +6x) v i2 +(3x 2 -4x+l)Ki]+(3x 2 -2x) K i2 } (1.3) 
dB 2 B?2-B?, 


-^r = (12x-6) Vii+(-12x+6)Vj 2 +(6x-4) K i( +(6x-2) K i2 (1.4) 

dB 2 


The values of Vj are known from the B-H curve and K; of (2) and (3) can be 
obtained by setting a constraint that the slope of K with respect to B 2 must to be continuous 
at the interval ends. For instance, if v-B 2 is divided into 6 sub-segments as shown in Fig. 
2, for two adjacent intervals from node i-1 to node i and from node i to node i+1, the 
gradient at the right end (x=l) is evaluated using (1.4) 
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dK 

-dB 2 Jj 


6Vj.i-6Vj+2Ki.i+4Ki 


(1.5) 


Similarly, calculating the quantity at node i at x=0 in the second interval 


dK 

-dB 2 Jj 


-6Vj+6v i+ |-4Ki-2K i+1 


( 1 . 6 ) 


Equating (5) and (6) 

Kj.i+4Kj+Ki + i = 3v i+ i-3Vj.i 


(1.7) 


For a v-B 2 curve with 6 subintervals, expanding (7) leads to 

K 1 + 4 K 2 +K 3 = 3 V 3 - 3 Vi 

k 2 +4k 3 +k 4 = 3v 4 -3v 2 
k 3 +4k 4 +k 5 = 3v 5 -3v 3 
K4+4K5+K6 = 3 v 6 - 3 v 4 
k 5 +4k 6 +k 7 = 3v 7 -3v5 


a set of 5 simultaneous linear equations with 7 unknowns. However, for the leftmost node 
at B=0, v can be easily obtained as a constant from B-H curve, leading to K]=0. For k 7 , 
different values must be tried until the series of Kj from (7) are positive and monotonic. 
With the above special treatment of leftmost and rightmost nodes, the unknowns are 
reduced to five 


4k 2 +K3 = 3V3-3V] 
k 2 +4k3+K4 = 3v 4 -3v 2 
K 3 + 41 Q+K 5 = 3v 5 -3v 3 
K1+4K5+K6 = 3 v 6 - 3 v 4 
K5+4K6 = 3V6-3v 4 -K 7 


Finally, use the calculated Kj in (1.2) and (1.3) to check the predicted the curve of 
H vs. B as illustrated in Figure B3. 



Silicon sheet steel 



Figure B3. Comparison of material sample test and curve modelling, 

circle material sample test 

square — curve modelling. 
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B.2 Iteration for Flux Distribution: Newton-Raphson Method 

To calculate the forces from a magnet accurately it is necessary to obtain solutions 
to the Poisson equation consistent with the magnetization characteristic of the metal. A 
Newton-Raphson iteration algorithm has been developed by Silvester et al.[B 1] and is 
suggested by many researchers working in this field due to its rapid convergence and 
unconditional stability. After evaluation of several iteration methods, the Newton-Raphson 
algorithm was found to converge faster and more stably than other methods, but it can not 
well constrain the saturated magnetic flux density because of the sharp increase of the 
reluctivity beyond the saturation point. However, with an undercorrecting parameter for 
the residual part of the iteration as well as an additional successive overrelaxation (SOR) 
form for a weighted combination of the reluctivity obtained at step n and step n+1, the 
Newton-Raphson algorithm offers a satisfactory solution. 

A summary of the iteration algorithm based on [Bl] is presented below. 

As presented in the previous report, the magnetic potential field will be obtained by 
minimizing the energy functional 

S(A> 

here W is the density of magnetically stored energy 

W: 

The requirement to minimize functional (2. 1) is equivalent to demanding 

^ n 

8A^ 0 ' < 2 -3) 

Expanding (2.3) will yield N simultaneous nonlinear equations whose solution A 
describes the desired result. 


To construct an iteration process, expand (2.3) in a Taylor series 
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(2.4) 


Neglecting those terms beyond the second derivative term will yield a matrix 


equation 


AA, 


JJSJ' (21) 

\3Aj5Aj /a+aa l^Aj / ( 


A+AA 


(2.5) 


Then a successive correction will be formed based on (2.5): 

A^A-Lf-S-r (2Ll 

\r)A|3Aj /a+aa l^Aj /a+aa 


( 2 . 6 ) 


To apply the iteration process (2.6) for finite element formulation, we calculate the 
first and second derivatives of (2.1) 


d§__ | / aw 
9Aj laAj 


_a 2 ^ 

3A; 


S_= ( d 2w 
idAj J dAjdAj 


j dn 


dQ 


(2.7) 


( 2 . 8 ) 


The field vector H is related to the flux density B by nonlinear reluctivity, 

H = v(B 2 ) B (2.9) 


Therefore (2.2) will become 

W=i-| v(B 2 ) d(B 2 ) 


4/ 


(2.10) 


Within an element, the potential is given by an interpolation form suitable for finite 
element calculations 


A= Yj AjNj 


( 2 . 11 ) 



Then the squared flux density in the element will be given by 


B 2 =22 A i A i VN i- VN j ( 2 - 12 ) 

' j 

Differentiating by the chain rule. 


— = v(B 2 ) X A k VNj- VN k (2.13) 

3Aj k 


3 2 W 

3Aj3Aj 


= vVNi-VNj + 2-d-X X A mAn(VN m -VNi)(VN n -VNj) 
dB 2 m n 


(2.14) 


Hence the residual vectors of the iterative process (2.6) will be formed by joining 
the individual element contributions formed according to (2.7) and (2.13) 




vVN k VNidQ- 


1 


JdQ 


(2.15) 




dQ= vVNj-VNj + 2 


I I A m A n | JMVN m VNi)(VN n VN,)d£i 
m n J R dB 2 


For a first-order triangular element, define 


I v VN VN j 2 1 VN m VNjd£2 = S mj 

Jr ‘ J /R 7r 


(2.17) 


and then the matrix contribution will be 
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|§S dn = '' S « + X^ISA m A„S imSjn 


The second term of (2. 19) can be simplified to 

X^X X A m A n S im S jn = UjUj 


where 


Ui = XSijAj 


So that the Jacobian matrix contributions are given, finally, by 

Pi ' = vSi ’ + W U ‘ u i 

Within every iteration, the reluctivity v is weighted by SOR: 


V (n + ')= V (n)-Ko( v *_ v (n)) 


(2.19) 


( 2 . 20 ) 

(2.21) 


(2.22) 


(2.23) 


where v* is calculated from the approximate magnetization characteristic. 
The iteration process is illustrated by the program flow chart of Figure B4. 
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B.3 Efficient Calculation of Force 

In static electromagnetic field analysis, the virtual work principle is commonly used 
to calculate magnetic force acting on a rigid, current-free and movable body: 


Fc = - 


aw 

as 


(3.1) 


where ds is an infinitesimal distance along the s-direction and W is the magnetic energy 
stored within the computation region. Therefore, to obtain magnetic force, the suspended 
body is subjected to a virtual displacement, then the field quantities are computed at both 
positions. A general procedure is first to solve the governing equation of a field problem 


with specified boundary conditions. For example, a static 2-D magnetic field will be 
expressed as a partial differential equation of Poisson's type: 


VvVA z (x,y) = J z . 


(3.2) 


Then the magnetic flux density will be obtained by calculating the curl of the vector 
magnetic potential 


B = V x A. ( 3 . 3 ) 

Magnetic energy finally is calculated by an integration over the volume containing 
the total distributed magnetic energy in terms of field quantities 

W 4|^ dv ' (3.4) 

In a numerical analysis for an electromagnetic field problem, the computation 
region may be bounded by irregular contours and may also involve several different 
materials. So a discretization of the domain may require many elements and the resulting 
force calculation would be very time-consuming. The problem becomes even more severe 
when an iterative method must be used to deal with nonlinearities such as the magnetization 
characteristics of different magnetic materials. To develop a more efficient and forthright 
algorithm for magnetic force calculation, Coulomb derived an elegant formula for 
implementing the virtual work principle without the need for a second field solution 
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[B3,B4], A mapping approach is used so that the calculations of differential terms of the 
field quantities have been switched to the calculations of differential changes of the 
coordinates in the Jacobian matrix. This approach not only saves the time that would be 
required for the second field solution, but may also improve the accuracy of the solution. 
It reduces the computer round-off error since the domain of the force calculation includes 
only the distorted elements during a virtual displacement. In the present research the 
principle of Coulomb's work is applied to derive the force calculation equations for two 
theoretical models: a linear solution for magnetic potential based on Laplace's equation, 
which was described in detail in a previous report, and a nonlinear Poisson's equation 
model. 

B.3.1. Force Calculation using Linear Potential Model 

With an assumption of infinite permeability of magnetic material, the computation region 
will only involve the air gap between the magnet poles and the suspended rotor. The 
fundamental equation of the source-free field then represents the curl-free nature of the 
magnetic field intensity 


V x H = 0.. 


(3.5) 


With the vector identity 


V x V(J) = 0,. (3.6) 

the field intensity can be written as the gradient of the magnetic potential 

H = -V<{>. (3.7) 

Since the divergence of magnetic flux density is zero everywhere 

V B =V p 0 H = 0. (3.8) 

substitution of (3.8) into (3.7) gives the Laplace equation 

V 2 <j> = 0. 


(3.9) 
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If the triangular linear interpolation formula 

<|>(x,y) = X Ni(x,y) 

i 

is used to approximate the potential field and then the energy functional 

i = ) iv N 2 dR, 

Jr 

is minimized as explained in the previous report, the solution of nodal potential values is 
found 

Jmin-fal. ®2, °nJ- 


(3.10) 


(3.11) 


( 3 . 12 ) 


Using the energy from equation (4), the force can be calculated from the virtual work 
principle of (1) , 

I ^dRl. (3.13) 

3s 3s t I 2(1,1 


where M is the total number of the elements. The Jacobian matrix can be used to transfer 
the global coordinates into local ones and to map a function between global and local 

coordinates 

j" f(x,y) dR= j f(x(Li,L 2 ),y(Li,L 2 ))|jldLidL 2 - (3.14) 
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The force then can be expressed as 


a function of flux density and local coordinates 


F s = 



[ Bi|jldL,dL 2 
2lK> 

Jn 



^2(10)' ' 2 ^to 9 s 


B: 




dLjdL 2 . 


J Q 


(3.15) 


The first term inside of the integral will be 


3(b 2 \_bT. 

asWo) 3s 

(3.16) 

With the equations of the mapping given as 


3<j)p 3<J>p 3x 3<|)p 3y 

3Li dx 3Li 3y 9Li 

(3.17) 

9d)R 3<t>p 3x . % dy 

, — + 5 

9L 2 dx 3L 2 dy dh 2 



or in a matrix form: 


' 5<f>p \ 


34>p \ 

9L, 1 

i = JI 

dx \ 

\ <% 

1 

\ 3<t>p 

\ 3L 2 j 


\ dy ) 


(3.18) 
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where the Jacobian matrix is 


dx dy 
j_j dLj 3Lj 
dx dy 

dh 2 dL 2 


(3.19) 


and (he derivatives lhC ,acob ™ 


matrix 


So from (1-14), 


d(j) 

H = -V<j> = .j-‘/ dL > 


dL? 


d(J) 

dH aj-i / a Ll 


a s | 


d<f> 

dh 2 


In order to eliminate the differential 


terms of <j> over L, the unit 


(3.20) 


(3.21) 


matrix is used 


ds 


ds 


ds 


= 0, 


(3.22) 


so that 


dJ 1 - , dJ 

— = -J‘ l -r' 

ds d s 


(3.23) 
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Substituting (1-17) into (1-16), we find 


m_ T _, 

9s 



(3.24) 


and then substituting (1-18) into (1-12), we finally have 





Bilfl-l 


apf 


dR . 


9s 9s 

Following a general routine of a calculation for the Jacobian matrix, we have 


(3.25) 


J = 


X 1 3 Y J3 
X 2 3 Y 23 


(3.26) 


where X,Y are nodal coordinates and the following forms are used to simplify the 
expression 


Xij = Xi-Xj, Yy = Yj-Yj. 


With an infinitesimal displacement of the nodal coordinates, we can write 

^2£li AYi3 

9J _ As As 

dS ^23 AY 23 

As As 


and 


9|j]_ A(X 13 Y 2 3-Yi 3 X 23 ) 
9s As 


(3.27) 


(3.28) 


(3.29) 


Geometrically, a suspension system undergoing a virtual displacement could be 
represented by three parts, a region of the magnet as a fixed part, a region of the levitated 
body as an entirely moving part and a free space region between the fixed and moving parts 
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as a distorted part. Equations (1-21) and (1-22) are discussed separately for each of these 
three different regions. 


A . For a fixed element 


AXjj = AYjj =0, 


(3.30) 


so 


dJ _ dfl n 

3s 3s ' 


(3.31) 


B . For an element within an entirely moving body, the element shape will not be 
changed at all, as illustrated in Fig.B5, 

(2) i 


(1) i 



Figure B5. An element in a rigid moving body is moved from location 1 to 2. 


AXy = X< 2) - xj 0 -(Xf> - xj !) = AX - AX = 0, (3.32) 

AYjj = Y^ 2) - Y^ -(Yj 2) - Yj !) = AY' - AY = 0, (3.33) 


where the superscript represents the moment before and after the displacement. For the 
same reason, 

A[XjjY n j-YijX nj ] = (X lj V B jP^X u Y^'M(Y| j X oj PMY i jX^'>) = 0, (3.34) 


since 


V (2)_ Y (1) V^-V* 1 * y(2)_ y (1) v(2)_y<0_0 
X ij -X ij ’ Y nj _Y nj’ X nj -A nj’ Y ij *0 


(3.35) 


So the same conclusion is reached: 
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^l-M = o (3.36) 

9s " 9s 

Applying the above results for (1-19), all those elements in the moving and fixed bodies 
can be neglected during the virtual displacement. Only the distorted elements need to be 
considered for the force calculation: 

Fs= I 

distor. 

J R 

Since the distortion of the free space is arbitrarily decided, we may just take one layer of 
the elements surrounding the moving part as shown in Fig.B6 where, for more accurate 
calculations, the diagonals of the gap 1 and 2 are symmetrically arranged. 



-B 1 J- 


— H + |j| _1 — 

9s 9s 


dR 


(3.37) 


Figure B6. Geometry of one layer of elements surrounding the rotor where, 

x fixed node 

o moving node. 
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If the first term of (1-23) is expanded to matrix form 


terml=| -B T J ' 1 — H 

ds 


A a =- 




[^x^y. 


^23 Y31 

X32 X 13 



,(3.38) 


where A a is the element area, all the quantities in term 1 are known except for the derivative 
of the Jacobian matrix. Since one rectangular shape consists of two triangular elements 
which have different ranges of fixed nodes and moving nodes, the matrix has to be 
determined separately according to the elements with odd or even numbers. An example of 
the nodal range is shown in Fig.B7. 


1 



Figure B7. The range of fixed and moving nodes of the elements. 


For element = 1,3,5, odd numbers, if s = x, since 



Ax Ax 

it is easy to determine that 

AX 13 _ q AX 23 _ AY 13 = AY 2 3 _ Q 
Ax Ax Ax Ax 


(3.39) 


(3.40) 
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So we have 


0 0 

dx ‘•" 10 - 


Repeating the above calculations, we obtain the matrix for s = y 


dJ = [ 0 0 

dx '■ 0 - 1-1 


(3.41) 


(3.42) 


The elements with even numbers and furthermore, the second term of (1-23), also can be 
handled with the same calculation above. The final forms of the force components along x 
and y-directions respectively as a summation of one layer elements are obtained as 


llr f^AxYsi+^yXnl+y^ylY 

1 (odd L 


13 


X " < J ) x(^xY23-M>yX32)+;^ < t ) x+ ( l>y)Y23 


J even 


3.43a 


F y =— ^|X <t>y( < t>xY 3 l+<|)yX 1 3)-^ 2 x+(l)y)x 1 3 

1 (odd L L J even 


"X < t > y( < l ) xY23+<l>yX32)+^(^x+<l>y)x23 


(3.43b) 


where D is the depth of the magnet. 

To examine the accuracy of the equations (1-24), the magnetic suspension force of an 
uniform magnet shown in Fig.B8 has been calculated and a comparison between analytical 
and numerical results is shown in Table Bl. 


MMF=1A x 400T 


y 


'T' 


0 



Figure B8. Test electromagnet 



B20 


Table Bl. Comparison of Analytical and Numerical results 



Force(N) 

y/c=0.0 

Force(N) 

yilc=0.7 

CPU 

Charge($) 

Analytical 

solution 

19.68087 

218.67681 



Numerical one 
layer solution 

19.68088 

218.67640 

1.80 sec. 

1.27 

Numerical whole 
region solution 

19.68109 

218.70340 

1.87 sec. 

1.30 


2. Force Calculation using Nonlinear Flux Model 

The methods described above will applied to the case where the flux is determined 
by Poisson's equation. The only difference is that the magnetic field is no longer current- 
free and the flux density will be the curl of the vector potential 


B = V x A. 


(3.44) 
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If there is only a z-component of current density, then 


B x 



0<f> 

dx 


(3.45) 


The computation region for the force calculation will be one layer of elements 
adjacent to the movable rotor and the final equations are expressed in the form of (1-27). 


F x =;p- {X [ < t ) y( < t ) x(^32+X 1 3)-<}>y( Y23+Y3 i)|-l-^{<t) x -4-<t)y)(Y23-Yi3) 

^ odd 


+ 


X [-<l) y (^yY23-<l)xX32)+^<l) 2 x+<t>?)Y23 


(3.46) 


F y= 2|io ^ ^ x ^ y ^ 23+ ^ 31 ^ x ^ 32+ X 13 ^^ x+ ^ y ^ 13 ~^ 23 ^ 

+ X [<!>x(<l)yY23-4>xX32)-^ 2 x+<t>y)x23] 


(3.47) 


A comparison of the results of this method to those of the original method, in 
which the energy contained in the entire solution region is used in calculating the virtual 
work involved, and the new one is displayed in Table B2. 


Table B2. Comparison of whole region and one layer calculations 



Force(N) 

y/c=0.0 

Force(N) 
yic= 0.7 
x/c=0.7 

CPU 

Charge($) 

One laye 
calculation 

19.59116 

186.24043 

42.99 sec. 

11.40 

Whole region 
calculation 

19.67500 

186.90760 

Imin38 sec 

26.79 
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Abstract 


Active magnetic journal bearings are increasingly being used in a wider variety of 
turbomachinery applications because their magnetic forces can be controlled and thus, 
can be used to minimize any rotor vibration. However, most magnetic bearing research 
has ignored the curved shape of the magnet and has consequently inappropriately 
modeled its forces. By using an accurate model of the two-dimensional forces of the 
curved magnets, this paper models the unbalanced rotating dynamic motion of the shaft 
which is being minimized by two opposed pairs of axis independent, proportional- 
derivative flux controlled magnetic bearings. The nonlinear dynamic behavior of the rotor 
is then examined using the resulting equations of motion via numerical simulation and the 
harmonic balance method. This dynamic analysis consists of examining the shaft's free 
vibration and its potential energy and maximum steady-state amplitude versus the natural 
frequency ratio for varying control and system parameters. 


1 
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Chapter 1 
Introduction 


A typical magnetic suspension system, see Figure 1.1, is comprised of two main 
parts, the actuator and the controller. The actuator, which often consists of two 
opposed pairs of magnets, see Figure 1.2, is responsible for levitating the shaft. While 
the controller, which consists of the feedback sensor(s), control electronics, and 
power amplifier(s), is used to regulate the forces of the actuator by using the signal 
from the sensor(s) and, depending on what is needed to minimize the shaft vibration, 
either increases or decreases the current supplied to each of the magnets. 

In recent years, this general type of active magnetic bearing has increasingly been 
used in a variety of applications. They are now found in centrifugal compressors, 
electric power plants, petroleum refining, machine tools, satellites and military 
weapons (O'Connor (1992)). This rise in popularity can mainly be attributed to two 
advantages that magnetic bearings have over the more conventional fluid film 
bearings. First, because in the magnetic actuator the rotor is levitated and no contact 
occurs between it and the magnets, there is practically no friction and thus, there is no 
wear and no need for lubricants. The other benefit of magnetic bearings is that they 
have the ability to minimize any shaft vibration. In recent years, this control of the 
shaft vibration has been the main topic of most of the magnetic suspension research 
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and has focused on two distinct subareas: 

1. The most suitable controller and control strategy 
2. Theoretically investigating the effectiveness of the control scheme. 

In trying to determine the optimal magnetic bearing control system, most 
researchers have experimented with digital and analog controllers and variations on 
proportional-derivative (pd) schemes. For example, Williams et al. (1991) 
experimentally compared the differences of using digital and analog pd controllers to 
reduce the amplitudes of a multimass flexible rotor. They found that in both control 
systems, increasing the stiffness tended to decrease the response at low speeds and 
increase the amplitudes at high speeds and increasing the damping minimized the 
maximum displacement for every forcing frequency. Williams et al. (1991) also 
concluded that including second derivative feedback in a digital pd controller gave 
the system a wider bandwidth while introducing integral feedback only affected the 
amplitude at the first critical speed. The possibility of using a digital pd controller 
was also experimentally examined by Scudiere et al. (1986); however, in their system 
the previous inputs were given uneven weight. The results showed that when the most 
recent inputs are heavily weighted the output responds slowly but is very smooth and 
that if the previous inputs are given less weight then the motion becomes less damped 
but responds much more quickly. 

The other popular magnetic bearing research topic is to theoretically investigate 
controlling the motion of a flexible rotor with a magnetic bearing. Lee and Kim 
(1992) designed a suboptimal output feedback controller based on a truncated modal 
equation of the distributed parameter system. This control scheme minimized the 
displacement at the first critical speed but was inadequate at other supercritical 
frequencies. Nonami et al. (1990) experienced similar problems at high frequencies as 
Lee and Kim (1992). Their flexible rotor model was based on the finite-element 
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method and the control scheme was formulated by using eigenvalue analysis. Maslen 
and Bielk (1992) also examined the stability of a rotor supported by a magnetic 
bearing. They took into account sensor-actuator noncollocation and controller 
bandwidth; however, the performance of the model was not tested. 

However, with the exception of Imlach et al. (1991) and Knight et al. (1992) there 
has been very little research completed which has measured the forces of the magnet. 
Imlach et al. (1991) compared predicted stiffness and force values in a closed-loop 
controlled magnetic bearing with their experimental counterparts. The agreement was 
good at small eccentricities; however, as the eccentricity increased the predicted and 
measured values began to diverge. Knight et al. (1992) also measured the forces from 
a magnet. They found that a single curved magnet produced forces in both the 
horizontal and vertical directions and that the horizontal force, assuming that the 
center of the magnet is situated directly on the y- or vertical axis (see magnet one in 
Figure 1.2), is proportional to the vertical force multiplied by the horizontal 
displacement. Most magnetic bearing papers, including Lee and Kim (1992), Nonami 
et al. (1990) and Maslen and Bielk (1992), have ignored this additional force 
component in their formulation of the vibration of the rotor. 

This thesis will show that including the often ignored normal force of a curved 
magnet in a magnetic bearing model will cause the motion of the unbalanced shaft to 
be significantly different than was previously estimated by earlier linear models. This 
will be accomplished by first using the conclusions of Knight et al. (1992) to derive 
brand new equations of motion for a shaft which is encircled by a flux controlled 
magnetic journal bearing, chapter 2. Then in chapter 3, the resulting coupled, 
nonlinear equations of motion and two techniques, numerical simulation (fourth-order 
Runge-Kutta) and the harmonic balance method, will be used to examine how the 
introduction of the normal force component changes the potential energy and 
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dynamic motion of the shaft for varying control as well as physical parameters. In the 
final section, chapter 4, the results from these two chapters will be briefly 
summarized and some suggestions will be made concerning future work. 
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Chapter 2 

Modeling Shaft Vibration 


The magnetic bearing model used in the upcoming analysis is comprised of a 
shaft and two opposed pairs of flux controlled magnets, see Figure 2.1. Each of the 



Figure 2. 1 : Magnetic bearing and shaft coordinate system. 


magnet pairs is identical and is subjected independently to proportional and 
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derivative axis control. In addition, the circular shaft is being externally rotated and 
its center of mass is not necessarily coincident with its actual geometric center. This 
rotating unbalance of the shaft, in addition to the forces of the magnetic actuator, 
cause the rotor to vibrate. The following section models these forces and the motion 
of the shaft, nondimensionalizes the resulting equations and analyzes appropriate 
parameter values. 

2.1 Dimensional Equation Derivation 

The magnets in Figure 2.1 are labeled one through four. In the initial portion of 
the equation derivation, the opposed pair two and four will be ignored and full 
attention will be given to the forces from the magnet pair of one and three. 


2.1.1 The Vertical Forces of Magnets One and Three 


The y-directional forces of these two magnets, F yl and F y3 , are approximated by 
one -dimensional magnetic circuit theory and are written as 


and 



aBjf 

Ho 



( 2 . 1 ) 


( 2 . 2 ) 


where B is the magnetic flux density of each magnet, a is the cross-sectional area of 
the pole face and (i 0 is the permeability of the free space. The total vertical force 

produced by the pair of one and three is simply the sum of (2.1) and (2.2), 


7 



C18 



Both flux densities in this expression are equivalent to the sum of their bias and 
control flux densities, denoted by subscripts b and c respectively, 

Bi=B lb +B lc (2.3) 

and 

B 3 = ®3b + ®3c» (2-4) 


and can be replaced to give 

(B lb + B, c ) 2 -(B 31) +B3 c ) 2 ]. (2.5) 

Since both magnets are identical, their bias flux densities must also be equivalent. 
Bib = ®3b = B b . The control flux densities of the two magnets however, will have 
opposite signs, B lc = -B 3c = B c , so that each magnet will force the shaft toward the 
origin. By replacing B lb , B lc , B 3b and B 3c in equation (2.5), with their equivalent, 
either B b or B c , expanding the bracketed terms and simplifying, F yl3 can be 

rewritten as 



Fy!3=— B b B c . (2-6) 

The control flux, B c , is a function of the vertical displacement, y, and velocity, y , the 
bias flux density, B b , the proportional control coefficient, k, and the derivative 
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control coefficient, y, and is equivalent to 


B c =-KB b y-yB b y. (2.7) 

Replacing B c in equation (2.6) with (2.7) yields the most useful representation of the 
vertical force caused by the magnet pair one and three, 

F yl 3 = — B,(icy + Ty). (2.8) 

M*o 

2.1.2 The Horizontal Forces of Magnets One and Three 

It was stated earlier, that Knight et al. (1992) found that a curved magnet 
produced a normal force (in the x-direction for magnets one and three and in the y- 
direction for magnets two and four) which was proportional to its normal 
displacement multiplied by the principal force (the force perpendicular to the normal 
force) of the same magnet; thus, the x-directional forces from magnets one and three 
can be written as 


F x! = a^Fyil 


and 

F x3 =a 3 x|F y3 


(2.9) 

( 2 . 10 ) 


where cq and cc 3 are the normal force proportionality constants and x denotes the 
horizontal component of displacement. The total x-directional force from these two 
magnets is simply obtained by summing (2.9) and (2.10), 
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F xl3= x (ai|Fyl| + a3| F y3|)- 


Because both magnets are identical, and therefore must have the same physical 
properties, oq and a 3 can be replaced by a to give 


F x13 



F x13 can be rewritten as 


F xi3 =—• ax(Bj+B|) 

M-o 

by substituting the absolute of equations (2.1) and (2.2) for F yl and F y3 . It is also 
useful to substitute (2.3) and (2.4) for Bj and B 3 , 



(B lb +B lc ) 2 + (B 3b 



( 2 . 11 ) 


This expression can be further reduced by replacing each flux density term with its 
respective equivalent, either B b or B c , and simplifying, 

F xi3=^-ax(BS+B?). (2.12) 

Finally, replacing the control flux density in (2.12) with (2.7) and expanding puts 
F x13 in its most practical form, 
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F x i3 =— axBj(l + K 2 y 2 +2icyyy + Y 2 y 2 ). 
M*o 


(2.13) 


2.1.3 All Other Forces Acting on the Shaft 

A similar procedure can be used to obtain the vertical and horizontal magnetic 
forces from magnets two and four. 


Fy24 = — ayB^(l + K 2 x 2 +2icyxx. + y 2 x 2 ) 
*^o 


and 


F X 24 = (kx + yx). 

M'o 


(2.14) 


(2.15) 


These two equations, along with (2.8) and (2.13), represent an accurate mathematical 
description of the forces from the magnetic actuator which act on the shaft. The shaft 
however, will also be subjected to an additional force which results from it being 
rotated and its mass being unevenly distributed. This additional load is analogous to 
rotating unbalance in rotor dynamics and its components can be written as 
meco 2 sincot, in the y-direction, and meco 2 coscot, in the x-direction, where m is the 
mass of the shaft, e is the distance from the geometrical center of the shaft to its 
center of gravity and co is the forcing frequency. By applying Newton’s second law of 
motion, £F = ma , in each direction, the nonlinear, coupled equations of motion 
which describe a shaft that is unbalanced, rotating and subjected to flux controlled 
magnetic fields are obtained, 


4a -^2 
my = B b 

Ho 


Ky + yy - — y(l + k 2 x 2 + 2 Tcyxx + y 2 x 2 ) 


+ meo) sin cot 


(2.17) 
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4a Ta2 
mx = B b 

Ho 


kx + yx - — x(l + K 2 y 2 + 2icyyy + y 2 y 2 ) 


-fmeco 2 coscot. (2.18) 


In both of these equations, note that all of the linear displacement and velocity 
terms are a direct result of the principal forces of the four magnets (the forces which 
can be modeled using one-dimensional magnetic circuit theory). In contrast, the 
normal forces are responsible for every coupling and nonlinear term in (2.17) and 
(2.18); thus, including the normal forces in the equations of motion of the shaft 
introduces the system to instabilities by decreasing its linear restoring forces, i.e. 
making it less stiff. 

2.2 Nondimensionalization of the Equations of Motion 


To make the upcoming analysis easier and the results more physically 
meaningful, (2.17) and (2.18) will be nondimensionalized by replacing the variables 
in both equations with the following, 


x = Xc 

II 

O 

t = T / co n 

co = n/co n 

K = K/c 

Y = r / cco n 

a = A/c 

e = Ec, 


where all capitals represent the nondimen sional quantities, c is the radial clearance 
between the rotor and the inner surface of the controlling magnets and the uncoupled 
natural frequency (a = 0) is equal to 
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CD 


n 



4aKB^ 


(2.19) 


Making these substitutions results in the following ordinary differential equations. 


cmcojY" 


ky + nr - — y(i + k 2 x 2 + iktxx ' + r 2 x /2 f 
^0 v 2 ' >) 

+ cmco 2 E£2 2 sin QT 


and 


cmco 2 X" = -Ap-^KX + TX' - y X(l + K 2 Y 2 +2KTYY' + T 2 Y' 2 ) j 


+ cmco 2 EQ 2 cosQT 


where a prime denotes differentiation with respect to nondimensional time. By 
dividing both of the two previous expressions by cmco^ and replacing the remaining 
squared uncoupled natural frequency terms in the resulting equations with (2.19) 
gives 


Y" = — 4am M- 0 B b r RY + py' _ A Y (i + k 2 X 2 + IKTXX' + T 2 X' 2 f 
4acmiqi 0 BbV 2 ' >) 


+ EQ 2 sinflT 


and 


X" 


— famH-oBb ( ^ + rx ,_^ x ( 1 + K 2 y 2 + 2KTYY' + T 2 Y' 2 )1 + EH 2 cosQT 
4acmK(i 0 B bV 2 ' ’) 


The most useful forms of the nondimensional equations of motion are obtained by 
canceling all like terms and replacing ck with K , 
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Y" = — 


■-UkY + TY' - ^ Y(l + K 2 X 2 + 2KTXX' + r 2 X' 2 )] + EH 2 sin QT (2.20) 


x" = 

Kv 


KX + EX' - yX(l + K 2 Y 2 + 2KTYY' + T 2 Y' 2 )j + E Q. 2 cosOT (2.21) 


Note that, in obtaining (2.20) and (2.21), c and the uncoupled natural frequency, 
co n , were chosen as the two variables which were used to aid in the 
nondimensionalization of equations (2.17) and (2.18). Both variables were selected 
over other parameters for different reasons. First, the choice of c allows the dependent 
variables in the governing equations of motion to become X, X , Y, Y and T. As a 
result, the numerical solution becomes more meaningful, because both X and Y, x / c 
and y / c respectively, are real physical properties and both are constrained to remain 
within the range of negative to positive one. In contrast, co n was chosen as a 
nondimensionalizing parameter because it will present the least amount of problems 
in the following analysis. Its only real drawback is that when the dimensionless 
normal force proportionality constant is varied, Q will not shift with it. Fortunately, 
because A is generally much smaller than K (their magnitudes will be discussed in 
the next section.), the actual shift in the forcing frequency ratio due to varying A is 
extremely small. The only other reasonable choice to use in place of co n is the linear 
natural frequency. 



Although this would allow the effect of A on £2 to be incorporated into the solution 
of (2.20) and (2.21), it would also mean that the equations of motion would have 
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1 / ^ J in P^ ace of both of the 1 / K and thus, some values would become 

cumbersome. For example, the uncoupled stiffness would change from a value of one 

,oK/ ( K -f} 

2.3 Appropriate Variable Ranges 

A significant part of the upcoming analysis is to examine how A , K , T and E 
affect dynamic behavior of the unbalanced rotating shaft. This will be accomplished 
by varying the four parameters within the following practical limits. 

The dimensionless normal force proportionality coefficient. A, for a single 
magnet was determined by Knight et al. (1992) to be approximately equal to 0.15. Its 
value for opposed magnets was less exact, so during the analysis, it will be varied 
from 0.05 to 0.25. The values of the dimensionless proportional and derivative 
control coefficients, unlike A , are less defined, although their selection does affect 
certain bearing properties and must be chosen accordingly. For example, magnetic 
bearings are often characterized by low damping; thus, to achieve an uncoupled 
damping ratio, r/2K, of 0.1 to 0.3 for a K equal to one, T must fall within the 
range of 0.2 to 0.6. K on the other hand is inversely proportional to the uncoupled 
damping ratio, so it would be advantageous to have large values of K which would 
decrease the damping in the bearing. Unfortunately, K is also equal to the inverse of 
the available displacement before the magnetic force becomes zero for T = 0, so to 
balance the two competing effects, K in the upcoming analysis will be limited to 
values under five. The fourth and final variable to consider, the nondimensionalized 
mass eccentricity, is dependent on the imperfection of the shaft, e, and the distance 
from the rotor to the magnet, c. Generally, e, the dimensional rotating mass 
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unbalance, is an extremely small number compared to c; however, c is solely 
dependent on the magnetic actuator design and in some instances may also be very 
small. It is unlikely, nevertheless, that in a magnetic actuator magnitude of e will 
approach that of c. A more reasonable situation is that the dimensionless mass 
unbalance is going to be equal to or less than 0.05; however, to be on the safe side, E 
will be given values as high as 0.15. 
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Chapter 3 

Results And Analysis 


3.1 Potential Energy 


If it exists (or can be derived), a useful tool in many dynamic analyses is a 
graphical representation of the potential energy, V, of the system. The dimensionless 
potential energy, which can be derived from the restoring forces in the conservative 
form of the governing equations, is 




k 2 x 2 y 2 

4 


(3.1) 


In addition to the vertical and horizontal displacements, note that V is a function of 
A and K . Consequently, fluctuations in either variable should lead to changes in the 
form of the potential energy. 

As A is increased from 0.05 to 0.25, with K held fixed, the potential energy 
wells, Figures 3.1a to 3.1c, become increasingly shallow. Thus, as expected, when the 
dimensionless normal force proportionality coefficient is increased the deviation of 
the potential energy from the linear case becomes more pronounced. This departure 
from the linear case is also increasingly evident at larger displacements. Also note 


17 



C28 





Figure 3.1: Potential energy wells where K=3 and a. A=0.05 b. A=0.15 c. A=0.25. 
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that in Figure 3.1 the potential functions, unlike some nonlinear systems, do not 

contain any unstable regions, i.e. a local maximum, and only have one equilibrium 

point, the origin, within the physical system. 

If K is now varied while A is held constant. Figure 3.2, the changes in the 

potential wells are similar to those observed in Figure 3.1. This may be somewhat 

surprising given the fact that typically, when something thought to be analogous to 

stiffness, such as the dimensionless proportional control coefficient, is increased, an 

increase not a decrease in stiffness should occur. However, it is important to 

remember that because of the nondimensionalization, the linear stiffnesses in (2.20) 

A 

and (2.21) are equal to 1 and thus, a reasonable increase in K will only affect 

2K 

the potential wells slightly. If K is now increased to twenty, Figure 3.3, there exists 
the possibility that the mass can enter an unstable area at high absolute values of X 
and Y. In the magnetic actuator, the likelihood that K will reach such a high 
magnitude is small; however, if the physical boundary of the dynamic system is 
increased, then the unstable equilibrium points corresponding to 


(X,Y) = (± 


2K- A 


AK 


2 »— 


2K- A 
AK 2 } 


will need to be avoided. 

3.2 Dynamic Motion of the Shaft 

In the following two sections, two different techniques will be used to exa min e 
the free and forced nonlinear dynamic behavior of the circular shaft and the effects of 
K , A , T and E on this motion. The initial analysis will be based on the Runge-Kutta 
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Figure 3.2: Potential energy wells where A=0.15 and a. K=1 b. K=3 c. K-5. 
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Figure 3.3: Potential energy well where K=20 and A=0.15. 

routine which is a highly accurate numerical method for solving differential equations 
and is described in Appendix A. In the second section, the results will be calculated 
using an approximate but sometimes extremely valuable analytical technique called 
the harmonic balance method. 


3.2.1 Runge-Kutta Method Free Vibration Results 

Often when investigating nonlinear systems, such as a pendulum undergoing 
large amplitude swings, it is useful to examine how the natural period changes with 
amplitude. Unfortunately, explicitly solving for the natural period in this case cannot 
be achieved because of the nonlinearity in the equations of motion of the shaft. 
However, by using numerical simulation, an accurate estimate of the natural period of 
the shaft as a function of the initial conditions can be obtained. An additional tool in 
the construction of these graphs is to use the symmetry of the system. For example, in 
Figure 3.3, the quartered section where both X and Y are positive, call it region one, 
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has the same shape as the other quarters of the well and within region one, its two 
halves are mirror images of each other. Because of this symmetry, only an eighth of 
the initial conditions need to be analyzed to gamer a complete understanding of the 
total behavior of the system. These initial conditions are defined as 


and 


Y 0 = Rsin({> 
X 0 = Rcos(J> 


where R is the nondimensional radial distance of the shaft and <{> is the angle 
measured from the x-axis to the radial distance vector, see Figure 3.4. 



Figure 3.4: Summary of natural period initial conditions. 

In Figure 3.5a, <J> is fixed at 22.5 degrees and the radius is increased from zero to 
one while both of the initial velocities are held at zero. For each increment in R the 
natural period ratio, x, for both displacements is determined by measuring the time 
between consecutive peak amplitudes and by dividing this result with the linearized 
natural period, 2k / ^l--^-,. Unlike the linear case, the x of both displacements is 

initial condition dependent and the values of the two natural period ratios diverge and 
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thus become more nonlinear as the radius increases. If the radius is now held fixed at 

0.5 and x for both X and Y is measured for each increment in <(), a graphical display 

of the angle versus the ratio is obtained, Figure 3.5b. Only when <)> is equal to 

— (2n -1) radians, where n is any integer, is there any agreement between the natural 
4 

period ratios of the two displacements. At every other angle the ratios are different 
suggesting that the undamped, unforced motion of the shaft will not repeat every 
period like the linear system. Instead, the behavior, see Figures 3.6a and 3.6b, will be 
analogous to beating, which occurs in harmonically forced linear systems where the 
natural and forcing frequencies are extremely close but not exact, and results in a 
trajectory which is not an ellipse or a circle but is closed and passes over a large 
portion of area. Figure 3.6c. 

3.2.2 Runge-Kutta Method Forced Vibration Results 

A useful tool for investigating how the nonlinearities affect the rotating dynamic 
behavior of the shaft over ranges of system parameters is to construct amplitude 
response curves for varying values of K , A , T and E. This graph plots the maximum 
steady-state nondimensionalized horizontal and vertical displacements versus the 
frequency ratio, Q and Figures 3.7 to 3.15 are obtained by numerically simulating the 
nondimensional equations of motion and by picking off both maximum 
displacements after 400 cycles (to allow for transients to decay). Unfortunately, 
contrary to their linear counterparts, the maximum amplitudes in nonlinear systems 
are initial condition dependent. Thus, to make the analysis easier, all of the following 
amplitude response curves will be constructed using the initial conditions 
X(0) = Y(0) = X(0) = Y(0) = 0, the rest state, unless stated otherwise. 

A comparison of Figures 3.7a to 3.7c, amplitude response curves where only K is 
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varied, yields some interesting results. First, as the dimensionless proportional control 
coefficient is increased the amplitudes grow, especially near resonance. This behavior 
may be somewhat surprising if K is again thought of in terms of stiffness or 
proportional feedback; however, because of the nondimen sionalization of (2.20) and 
(2.21), K is actually more closely related to the inverse of the uncoupled damping 
ratio, T/2K, in each equation; therefore, the growing amplitudes near resonance 
make sense for increasing values of K. The most surprising feature of graphs 3.7a to 
3.7c is the difference in the maximum amplitudes of X and Y for K = 5. In the 
analogous uncoupled linear case, A = 0, the maximum amplitudes are identical, 
Figure 3.8; thus, the amplitude split in Figure 3.7c can be attributed to the 
nonlinearity of the equations of motion. It is not completely surprising that this 
behavior only occurred when K was large. For increased values of K, larger 
amplitudes are expected and earlier, while examining the potential wells, larger 
displacements meant the nonlinearity of the equations became more pronounced. 

The time series, trajectories and phase projections of the nonlinear and linear 
K = 5 cases at comparable frequency ratios, Figures 3.9 and 3.10 respectively, yield 
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Figure 3.9: Numerical steady-state a. Time series b. Trajectory c. Phase portrait 
where K=5, A=0.15, T=QA, E=0.1 and 0=0.9925. ‘ 
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Figure 3.10: Numerical uncoupled steady-state a. Time series b. Trajectory c. Phase 
portrait where K=5, A=0, r=0.4, E=0.1 and £2=1. 
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some further interesting behavior. (Note that Figures 3.9 and 3.10 are at different 
values of , this is necessary because a was not included in the selection of co n .) As 
expected, in the linear time series, Figure 3.10a, the maximum amplitudes are 7C / 2 
out of phase. This however is not true in the nonlinear case. Figure 3.9a.; therefore, 
when one component of displacement is at its maximum the other has a nonzero 
value and the path resembles an angled ellipse. In addition, the maximum radial 
amplitude in the nonlinear case, Figure 3.9b, is larger thus increasing the likelihood 
of the rotor striking the magnetic actuator. It should be noted that, although the actual 
maximum amplitudes and the transient motion in Figures 3.7 to 3.10 and some of the 
subsequent graphs are well outside the physical limitations of the bearing, it is 
beneficial to examine these cases because they provide a better dynamic 
understanding of the system and may prove even more useful if the amplitude 
restriction is somehow changed. 

In Figures 3.11a to 3.11c, every variable except the normal force proportionality 

constant is fixed. In each of the three graphs the amplitudes tend to decrease only 

slightly and the peaks and the entire plots tend to be slighdy shifted to the left for 

increasing A . Once again this shift in amplitude is purely a function of what natural 

a 

frequency is used in the nondimensionalization of the equations of motion. If k — — 

replaced K in equation (2.19) then no shifts would occur in Figures 3.11a to 3.11c. 
The major difference between the three graphs is, in the third graph. Figure 3.11c, 
another amplitude split occurs. However, unlike the two distinct maximum 
amplitudes in Figure 3.7c, these different amplitudes result because for greater values 
of A the system will deviate more from the linear case, recall Figures 3.1a to 3.1c. 

Unlike K and A, T must be decreased for the maximum X and Y displacements 
to split. Actually, the effect of the dimensionless derivative control coefficient on the 
amplitude response curves, Figure 3.12a to 3.12c, is extremely similar to the 


30 


















C43 


influence of 1/ K on Figures 3.7a to 3.7c; in that any increase in T results in smaller 
maximum amplitudes especially near resonance. This relationship between the two 
variables is not surprising because the uncoupled damping ratios in the equations of 
motion are equal to T / 2K . In Figure 3.12a, unlike any other plot up until this point, 
an additional solution occurs directly to the left of the amplitude split. This tiny 
branch of solutions was found using different initial conditions from the traditional all 
zeros. Instead, the frequency ratio was slowly decreased from 2.5 and the starting 
values for each decrease in were obtained using steady-state data from the 
preceding frequency ratio. Using this method, the alternate solution was tracked until 
approximately Q. = 0.8875. At that point, the only solution that could be located was 
the one where the two maximum amplitudes agree. It is also likely that an unstable 
curve exists from the lowest frequency ratio of the alternate solution to the initial 
point of the amplitude split, Q. = 0.9125. Tracking this unstable motion using present 
path-following algorithms is however beyond the scope of the current work. Also, 
because of the high complexity of investigating every initial condition, no specialized 
attempt will be made to locate other possible stable solutions. 

A more complete picture of how A, T and K, because of its close inverse 
relationship to T, change the amplitude response curves and thus the dynamic motion 
of the rotor, is shown in Figure 3.13. The results confirm the previous analysis. An 
increase in the normal force proportionality constant will cause a slight decrease in 
amplitude but an increased probability that the maximum amplitudes will not be 
identical near resonance and that an additional solution will appear before resonance. 
Decreasing T as well as increasing K mainly cause an increase in the amplitude near 
resonance and an increase in the likelihood of an amplitude split. 

From a design standpoint it is also important to investigate how the eccentricity of 
the shaft affects the amplitudes at given frequency ratios. Figures 3.14a to 3.14c. 
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Figure 3.13: Numerical amplitude response curves for K=3, E=0.1 and varying 
A and F where a dashed line represents X, a solid line represents Y, and a single 
solid line represents both X and Y. 


34 
















C46 



Figure 3.15: Numerical amplitude response curve where K=3, A=0.175, T=0.2 
and E=0.2. 


Because E is proportional to the external rotating unbalance forcing amplitude, any 
increase in it results in an almost proportional increase in maximum displacement for 
a given £2. Thus, as E gets larger, the displacements increase and an amplitude split 
becomes more likely, Figure 3.14c. If the initial forcing phase of the system is now 
varied, the results in Figure 3.14 and all of the other previous graphs will remain 
generally the same. The one exception is that for some initial forcing phases the 
amplitudes and phase shifts of X and Y may be interchanged. 

Other interesting, nonlinear behavior besides the amplitude jump can occur with 
the governing equations. For instance, in Figure 3.15, an amplitude response curve 
where K = 3, A = 0.175, T = 0.2 and E = 0.2, at approximately £2 = 0.9 and again 
at £2 = 0.9875 there exists multiple maximum amplitudes for the given frequency 
ratios. At both frequency ratios, the maximum X and Y displacements fluctuate from 
period to period, Figures 3.16a and 3.17a, suggesting either chaotic or quasi-periodic 
motion (Thompson and Stewart (1986)). In order to correctly classify this behavior it 
is necessary to complete a spectral analysis of the time series data using the Fast 
Fourier Transform (FFT). The FFT converts 2 N , where N is any positive integer, 
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Figure 3.16: a. Numerical time series b. FFT of X c. FFT of Y for K=3, A=0.175, 
r=0.2, E=0.2 and £2=0.9. 
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Figure 3.17: a. Numerical time series b. FFT of X c. FFT of Y for K=3, A=0.175, 
r=0.2, E=0.2 and 0=0.9875. 
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equally spaced discrete data points, making sure that there are at least two points per 
highest forcing period, to a frequency based domain using the Fourier Transform. 
After the data has been changed to its new domain, it will only have significant 
amplitudes at various combinations of its forcing frequencies. Thus, with the FFT, it 
is possible to classify any motion as periodic, one forcing frequency, chaotic, infinite 
forcing frequencies (also known as a broad-banded solution), or, as in the case of 
Figure 3.16 and 3.17, quasi-periodic, two or more driving frequencies. 


3.2.3 Harmonic Balance Method Forced Vibration Results 

The dynamic behavior of the circular shaft can also be analyzed using an 
analytical method called harmonic balance. This technique predicts the steady-state 
response of the system by assuming solutions of the form 

Y = CcosfiT + Dsin£2T (3.4) 

and 

X = GcosfTT + HsinQT (3.5) 


where C, D, G and H are unknown constants. Substituting (3.4) and (3.5) and their 
first and second derivatives with respect to nondimensional time, T, into (2.20) and 
(2.21) yields the first two equations shown in Appendix B. Both equations can be 
simplified by expanding, replacing the higher powered terms with their appropriate 
trigonometric relations, 


cos 3 0 = — (3 cos 0 + cos 30) 


sin 3 0 = i(3sin0 + sin 30) 
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cos 2 9 sin 9 = (sin 9 - sin 39) cos9sin 2 9 = (cos 9 - cos 3 9) 

and neglecting the higher harmonics. See Equations B.3 and B.4 in Appendix B. 
Next, the necessary equations needed to solve for the four unknowns, C, D, G and H, 
result by separately grouping all of the cosfit and sin£2t terms in B.3 and B.4 and 
setting them equal to zero (balancing the harmonics). 

However, due to the coupling and higher-order terms in the resulting gove rnin g 
equations, B.5 to B.8, the constants cannot be solved for analytically. Instead, they 
must be determined by using a numerical procedure called Newton's method for 
nonlinear systems. The method of solution for this technique is analogous to the more 
well known Newton-Raphson method. Both methods use an initial guess of the 
unknown(s), the partial derivatives of the function(s) and iteration to converge 
quadrically on a local solution. In cases where multiple solutions exist, however, the 
result can be solely dependent on the starting values of the unknowns and thus, 
different initial guesses may yield different results. This however makes it possible to 
locate a variety of approximate analytical solutions for given K , A , T and E by 
simply choosing different initial guesses for C, D, G and H. 

Since the approximate analytical approach ignores all higher harmonics, it is 
necessary to examine if, in this case, the technique yields valid results. This can be 
accomplished by simply comparing graphs which are constructed using the 
approximate analytical approach, Figures 3.18a and 3.19a, with graphs which are 
obtained using numerical simulation, Figures 3.18b and 3.19b. The first two graphs 
plot the displacement versus the frequency ratio, Q, at any steady-state 
nondimensional time n7l, where n=0,2,4,..., and Figures 3.19a and 3.19b plot the 
velocity versus the Cl at the same time, n7t. Both sets of graphs are almost identical, 
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Figure 3.20: Harmonic balance amplitude response curve where K=3, A=0.15, T=0.2 and E=0.1. 


suggesting that the approximate analytical approach does provide an accurate 
description of the steady-state dynamic behavior of the shaft. As an additional check, 
an amplitude response curve with the same parameters as Figure 3.12a can be 
constructed using the approximate analytical approach, Figure 3.20. A quick 
comparison of these graphs continues to strengthen the validity of the approximate 
analytical approach because the plots are almost identical with the exception of a 
third solution near resonance in Figure 3.20. The absence of this additional solution 
in Figure 3.12a can be simply attributed to the fact that it is probably unstable and 
numerical simulation can only be used to locate stable solutions. 

3.2.4 Limiting Shaft Displacement for Varying Parameters 

Ideally, the unbalanced shaft displacement should be minimal and at worst it 
should never come into contact with the bearing. Thus, for certain values of K , A 
and T it is beneficial to know what nondimensionalized eccentricities will cause a 
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specified unwanted displacement and what effect varying each parameter has on these 
eccentricity values. Figure 3.21a, which was constructed using the approximate 
analytical approach, is a plot of the lowest eccentricity value which causes the 
amplitude of either X or Y, X max , to exceed 0.4, 0.7 and 1.0 versus the frequency 

ratio for K = 3, A = 0.15 and T = 0.4. As expected, the allowable nondimensional 
eccentricity for each specified amplitude is smallest near resonance and largest as 
£2 = 0 is approached. Somewhat surprising however is that below £2 h 0.75, the 
curves are identical. At these small frequency ratios, the lowest eccentricity value 
corresponds to the smallest E which causes the amplitudes to split. In each case 
below £2 = 0.75, the amplitude of X or Y for these eccentricities is always greater 
than one and as a result, the curves coalesce. The final noteworthy point of Figure 
3.21a is the presence of a kink, which is denoted by an arrow, in the X max = 1.0 case. 

This bend only occurs in large maximum amplitude cases because for large 
displacements, the amplitude response curves will not be smooth, as indicated by 
Figures 3.7c, 3.12a and 3.14c. 

The next three graphs, Figures 3.21b, 3.22a and 3.22b, continue to have the 
lowest E and £2 as the abscissa and ordinate variables; however, the amplitude limit 
is no longer the third variable. Instead, either K , A or T , depending on the plot, 
becomes the last parameter. As for the maximum allowable amplitude, it has been set 
at 0.4 in each graph because, typically, the displacements should be kept as small as 
possible. In the first of the three figures, the nondimensional proportional control 
coefficient has been increased from one to five. As anticipated, when £2 = 1, the 
smaller the K , the larger the eccentricity needs to be for either of the displacements 
to become 0.4, recall Figures 3.7a to 3.7c. It also makes sense that at high values of 
£2 there is very little difference in the required E; however, according to the 
amplitude response curves this same behavior is also expected at low frequency 
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Figure 3.23: Harmonic balance amplitude response curves where A=0.15, r=0.4, 
E=0.1 and K is varying. 


ratios. But as stated in the previous paragraph, at small frequency ratios the lowest 
value of the nondimension al eccentricity is dependent on the lowest E at which the 
amplitude split occurs. And for varying K, these values are different, Figure 3.23. In 
examining Figure 3.23, it is also important to recognize that for increasing values of 
K greater than five, there exists a range of frequency ratios where the amplitudes 
decrease and a larger E would be required to cause a specified displacement. 
Although this only occurs for values of K which induce a split, the entire behavior is 
completely contrary to any other previous results and would actually mean that if 
K = 6 and K = 7 were included in Figure 3.21b their solutions would cross. In Figure 
3.22a, the derivative control coefficient has replaced K as the third parameter. 
Because the effect of T on displacement is almost identical to the influence of 1 / K 
on the same variable, see Figures 3.7 to 3.12, Figures 3.21b and 3.22a exhibit similar 
results for the same reasons. 

Finally, the effect of the normal force proportionality constant on the lowest 
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eccentricity value which causes the maximum displacement of either X or Y to 
exceed 0.4 is examined in Figure 3.22b. For most of the frequency ratios, any change 
in A results in little or no change in the lowest eccentricity value. This agrees with 
the results from the earlier amplitude response curves, Figures 3.11a to 3.11c, where 
an increase in A meant little change in amplitude. However, once again, at small 
frequency ratios the lowest value of eccentricity does not relate to the amplitude 
response curves because, as in the case of K and I\ the value at which the split 
occurs varies for a changing A . 

Each of the graphs, 3.21a, 3.21b, 3.22a and 3.22b, were constructed using the 

approximate analytical approach. Although it has been shown that this technique is 

sound when examining steady-state behavior, it does ignore transient motion which 

may be significant especially if the rotor is started from rest. Thus, if escape from the 

2K- A 


potential well, defined as the magnitude of X and Y exceeding 


AK" 


at the same 


nondimensional time, T, occurs before the onset of steady-state behavior, then 
Figures 3.21 to 3.22 are rendered inappropriate. It should be noted that the transient 
behavior will actually cause the shaft to strike the bearing before escape can occur, 
unless K is unreasonably high, thus changing the equations of motion of the shaft, 
see Figures 3.2 to 3.3. This modification of the equations has been ignored in the 
upcoming graphical analysis; however, for a system with larger physical bounds the 
forthcoming figures are appropriate. 

Figures 3.24 and 3.25 display the nondimensional eccentricities, for a range of 
frequency ratios, which cause escape during the first one hundred dyn ami c cycles, 
after which transient motion has been assumed to have died out. (This is usually but 
not necessarily the case.) Comparing the results obtained in Figure 3.22a with Figures 
3.24a and 3.25, it is evident that the lowest eccentricity values which cause escape 
during transient motion are always larger than the values which cause the steady-state 
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amplitude to exceed 0.4. Thus, Figures 3.21 to 3.22 are valid. 

The Figures 3.24a and 3.25 are themselves noteworthy. First, in both graphs there 
is not a solid line which separates escape eccentricities from nonescape eccentricities 
and the boundaries which do exist are fractal (Crilly et al. (1991) and Feder (1988)). 
For instance, in Figure 3.24 as the resolution of the graph increases the boundaries 
change and if these escape boundaries are zoomed in upon further and further then 
the fractal behavior of the plot will continue infinitely. (It is possible to calculate the 
noninteger dimension of these boundaries, however, due to the time constraints of 
this thesis it was not attempted.) The actual shape of either graph is extremely 
complex as well, with gray areas, denoting escape between four cycles and one 
hundred cycles, often occupying regions which are completely surrounded by black 
sections, escape within the first three cycles, see Figure 3.24b. Yet each plot tends to 
exhibit expected behavior. For example, when the derivative control coefficient is 
increased the shape of either graph generally remains the same except the escape 
eccentricity decreases for all frequency ratios. Also, the minimum value of E in both 
graphs occurs near resonance, and at low frequency ratios, extremely high 
eccentricities are necessary for escape. Varying K and A does not produce any 
unexpected results either, Figure 3.26. Once again, the graph tends mainly to shift 
downward. The one unexpected and definitely nonlinear quality of Figures 3.24 to 
3.26, is the fact that for increasing T, thus damping, the minimum E tends to move 
towards a lower frequency ratio. This agrees with the nonlinear large amplitude 
motion of Figure 3.23; however, in linear rotating unbalance and even in the 
magnetically controlled shaft where the damping ratio is sufficiently high enough that 
the system is close to being linear, Figure 3.27, the maximum steady-state amplitude 
tends to move to the right for increasing T. 
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Chapter 4 

Conclusions and Further Work 


Recent active magnetic bearing research determined that curved magnets produce 
forces which are two-dimensional. Using this conclusion, this dissertation has 
accurately modelled the proportional-derivative flux controlled forces of the magnetic 
bearings and the motion of an unbalanced rotating shaft which is being minimized by 
these magnetic forces. The resulting governing equations are coupled and nonlinear 
and exhibit important atypical dynamic behavior such as quasi-periodic motion and 
hysteresis. In addition, the nonlinear terms also introduce the possibility of having 
unstable motion within the actuator and it is also possible that near resonance and 
with either large displacements or a large normal force proportionality constant, that 
the steady-state maximum amplitudes may not be identical and multiple different 
larger amplitude solutions may exist for a given set of parameters. Variations in any 
one of these parameters tends to affect the dynamic motion differently; however, 
generally, any increase in the nondimensional proportional control coefficient, K , or 
the nondimensional rotating unbalance, E, or decrease in the nondimensional 
derivative control coefficient, T, will cause an increase in amplitude and thus, 
increase the possibility of introducing the unwanted nonlinear behavior. An increase 
in A , the dimensionless normal force proportionality constant, will not increase the 
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amplitude; however, it will make the nonlinear larger amplitude motion more 
probable. 

Further magnetic bearing research still needs to be completed in two areas. First, 
because an increased nondimensional normal force proportionality constant can 
introduce unwanted nonlinear behavior without any change in amplitude, it is 
necessary to further examine its specific value in opposed pairs of magnets. In 
addition, the previously presented results need to be experimentally verified. 
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Appendix A 

Runge Kutta Method 


The fourth-order Runge-Kutta method is the most popular numerical technique 
that is used to solve nonlinear ordinary differential equations because it strikes a 
balance between the computational complexity and accuracy of the solution. The 
derivation of this method can be found in any numerical analysis text and yields the 
following formulas for two coupled, first-order differential equations, v = f(t,v,x) 
and x = g(t,v,x). 


Vn + i=v n +^(K,+2K 2 +2K 3 +K 4 ) 

(A.l) 

x n+l = x n + ^( L l + 2L 2 + 2L 3 + L 4) 

(A.2) 


where h is the time step, n is any integer, v 0 and x 0 are known and 


^1 f(tn» v n>*n) 

__ Y h K, LA 
K 2 = f[t n+ -,v n+ -^,x n+T J 


K, =f 


h K 2 L 2 

0 » x n ^ 
V 2 2 2 
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K 4 =f(t n + h,v n +K 3 ,x n +L 3 ) 

4 = g(t n >v n ,x n ) 

r ( h K, L,^ 

2 g[t n + 2 ,v n + 2 ,x n + 2 J 

T r h K 2 L 2 "| 

3 g[t„ + 2 ,v n + 2 ,x n + 2 J 

L 4 =g(t n + h.v n +K 3 ,x n + L 3 ) 


When the two coupled equations of motion are second-order, it is necessary to put 
each into state space form, 


v = w 

w = f(t,v,w,x,y) 


x = y 

y = g(t,v,w,x,y) 


The same general theory can now be used to solve these four coupled first-order 
equations, 


V „1 = v„ + |(W, + 2 W 2 + 2W 3 + W„) 

(A.3) 

w„ + i = w„ + V, + 2K,+2K, + K 4 ) 
0 

(A.4) 

^l=x„+|(Y 1 +2Y 2 +2Y 3 +Y 1 ) 

(A.5) 

yn+1 = Yn + ^( L l + 2L 2 + 2L 3 + l 4 ) 

(A.6) 


except in this case w 0 and y 0 are also known and 
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W 1= w n 

Ki =f(t n ,v n ,W 1 ,x n ,Y 1 ) 

Y i = y n 

Lj — g(t n , v n , W2,x n , Yj) 

W 2 =w i + |k 1 

L 2 = g(t„ + ^,V„ + 1 W„ W 2 ,x„ + 1 Yl Y 2 j 
Y 2 = y»+^L! 

K2 = f(t„ + |,v n +|w 1 /W 2 ,x 1 ,+|Y I ,Y 2 ) 
W 3 =w„+|k 2 

K 3 =f[t„+|,v„+|w 2 ,W 3 ,x„ + iY 2 ,Y 3 ) 

xx h T 

y 3 = y n + 2 L 2 
W 4 =w n + hK 3 

K 4 = f(t n + h,v n + hW 3 ,W 4 ,x n + hY 3 ,Y 4 ) 
Y 4 =y n + hL 3 

L 4 = g(t n + h,v n + hW 3 ,W 4 ,x n + hY 3 ,Y 4 ) 
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Appendix B 

Harmonic Balance Method Equations 


-CG 2 cosGT -DC! 2 sinGT + r(-CGsinGT + DG cosGT) +K(C cosGT + Dsin GT) 

- y (C cos GT + D sin GT)(1 + K 2 (G cos GT + H sin GT) 2 

+2KT(G cos GT + H sin GT)(-GG sin GT + HQ cos GT) + T 2 (-GG sin GT + HG cos GT) 2 ) 

= EC! 2 sinGT (B.l) 


— GG 2 cos GT - HG 2 sin GT + T (-GG sin GT + HG cos GT) + K(G cos GT + Hsin GT) 

(G cos GT + Hsin GT)(1 + K 2 (C cosGT + D sin GT) 2 

+2KT(CcosGT + DsinGT)(-CGsinGT + DGcosGT) + r 2 (-CGsinGT + DGcosGT) 2 ) 

= EG 2 cosGT (B.2) 


— CG 2 cos GT - DG 2 sin GT + r(-CG sin GT + DG cos GT) + K(C cos GT + D sin GT) 
-y((CcosGT + DsinGT) + K 2 (^CG 2 +^-CH 2 +iDGHjcosGT 

+Q-DG 2 +|DH 2 +^CGHjsinGT) + 2KrG(^-^DG 2 +^DH 2 + ^-CGh| cosGT 

+(- - CG 2 + - CH 2 - - DGh] sin GT) - T 2 G 2 (f- CG 2 + - DG 2 - - DGH 
{ 4 4 2 J \4 4 2 


+|-CH 2 +-DH 2 --CGH 
.4 4 2 


>= 


EG 2 sinGT 


sinGT 

(B.3) 
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-GQ 2 cos QT - HQ 2 sin QT + T(-GQ sin QT + HQ cos QT) + K(G cos QT + H sin QT) 

— 4 ((G cos QT + H sin QT) + K 2 (f— C 2 G + — D 2 G + — CDh) cos QT 
2 V4 4 2 J 

+^C 2 H + |D 2 H + ^CDGjsinQT) + 2KTQ(^C 2 H + iD 2 H + ^CDGjcosQT 

+ [-i c2 G + ±D 2 G - IcDHjsinQT) - r 2 a 2 (jV G + | G 2 H - loDHjsinQT 

(B.4) 


+|4-D 2 G + -D 2 H--CDg 1)) = EQ 2 cosQT 
V4 4 2 J 


(B.5) 

-CKQ 2 + CK + Dra-^C-y DGHK 2 - CG 2 K 2 - — CH 2 K 2 - — CGHKXQ 

2 4 8 8 2 

+^DG 2 KrQ-^DH 2 KTQ + ^DGHr 2 Q 2 -yCG 2 r 2 Q 2 -^CH 2 r 2 Q 2 = 0 


(B.6) 

-DKQ 2 + DK - Cm - -^-CGHK 2 - ^DG 2 K 2 - — DH 2 K 2 + — DGHKTQ 

2 4 8 8 2 

+~CG 2 KTQ - 4cH 2 KTQ + -CGH T 2 Q 2 - — DG 2 r 2 Q 2 - -DH 2 r 2 Q 2 = EKQ 2 
4 4 4 8 8 

(B.7) 

-GKQ 2 + GK + HTQ -^G- — CDHK 2 - — C 2 GK 2 D 2 GK 2 - — CDGKTQ 

2 4 8 8 2 

+4c 2 HKm--D 2 HKm + -CDHT 2 Q 2 --C 2 Gr 2 Q 2 D 2 Gr 2 Q 2 = EKQ 2 
4 4 4 8 8 


(B.8) 

-HKQ 2 + HK- Gm - — H - — CDGK 2 - — C 2 HK 2 - — D 2 HK 2 + — CDHKTQ 

2 4 8 8 2 

+4c 2 GKm- — D 2 GKTQ + A C DGr 2 Q 2 c 2 ht 2 q 2 --d 2 ht 2 q 2 = 0 

4 4 4 8 8 
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